/* nag_estimate_agarchI (g13fac) Example Program.
 *
 * Copyright 2000 Numerical Algorithms Group.
 *
 * Mark 6, 2000.
 *
 *
 * Mark 8 revised, 2004
 *
 */

#include <nag.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <nagg05.h>
#include <nagg13.h>

int main(void)
{

  Integer exit_status=0, i, ip, iq, isym, j, k, maxit, mn, npar, nreg, nt;
  Integer num, num_startup, seed, tdc, tdx;
  NagError fail;
  Nag_Garch_Est_Initial_Type est_opt;
  Nag_Garch_Fcall_Type fcall;
  Nag_Garch_Stationary_Type stat_opt;
  double *bx=0, *covar=0, *cvar=0, *et=0, fac1, gamma, hp, *ht=0, lgf, mean;
  double *param=0, *rvec=0, *sc=0, *se=0, *theta=0, tol, *x=0, xterm, *yt=0;

  INIT_FAIL(fail);
  isym = 1;
  gamma = -0.3;
  nreg = 2;
  ip  = 1;
  iq = 1;
  mn = 1;

  npar = iq + ip + 1;
  num = 1000;
  nt = 3;

  tdc = npar+mn+isym+nreg;
  tdx = nreg;

#define YT(I) yt[(I)-1]
#define THETA(I) theta[(I)-1]
#define SE(I) se[(I)-1]
#define SC(I) sc[(I)-1]
#define RVEC(I) rvec[(I)-1]
#define PARAM(I) param[(I)-1]
#define HT(I) ht[(I)-1]
#define ET(I) et[(I)-1]
#define BX(I) bx[(I)-1]
#define CVAR(I) cvar[(I)-1]
#define X(I,J) x[((I)-1) * tdx + ((J)-1)]
#define COVAR(I,J) covar[((I)-1) * tdc + ((J)-1)]



  if ( !(bx = NAG_ALLOC (nreg, double))
       || !(covar = NAG_ALLOC ((npar+mn+isym+nreg)*(npar+mn+isym+nreg), double))
       || !(et = NAG_ALLOC (num, double))
       || !(ht = NAG_ALLOC (num, double))
       || !(param = NAG_ALLOC (npar+mn+isym+nreg, double))
       || !(rvec = NAG_ALLOC (40, double))
       || !(sc = NAG_ALLOC (npar+mn+isym+nreg, double))
       || !(se = NAG_ALLOC (npar+mn+isym+nreg, double))
       || !(theta = NAG_ALLOC (npar+mn+isym+nreg, double))
       || !(x = NAG_ALLOC (num * nreg, double))
       || !(cvar = NAG_ALLOC (nt, double))
       || !(yt = NAG_ALLOC (num, double)))
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }


  Vprintf ("nag_estimate_agarchI(g13fac) Example Program Results \n\n");
  seed = 11;

  mean = 3.0;

  if (nreg > 0)
    {
      for (i = 1; i <= num; ++i)
        {
          fac1 = (double) i *0.01;
          X (i, 1) = sin (fac1) * 0.7 + 0.01;
          X (i, 2) = fac1 * 0.1 + 0.5;
        }
      BX (1) = 1.5;
      BX (2) = 2.5;
    }


  PARAM (1) = 0.15;
  PARAM (2) = 0.1;
  PARAM (3) = 0.8;
  PARAM (4) = 0.1;

  fcall = Nag_Garch_Fcall_True;
  /* nag_random_init_repeatable (g05cbc).
   * Initialize random number generating functions to give
   * repeatable sequence
   */
  nag_random_init_repeatable(seed);
  num_startup = num;
  /* nag_generate_agarchI(g05hkc).
   * Univariate time series, generate n terms of either a
   * symmetric GARCH process or a GARCH process with asymmetry
   * of the form (epsilon_(t~-~1)~+~gamma~)^2
   */
  nag_generate_agarchI (num_startup, ip, iq, &PARAM(1), gamma, &HT(1), &YT(1),
                        fcall, &RVEC(1), &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_generate_agarchI (g05hkc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  fcall = Nag_Garch_Fcall_False;
  /* nag_generate_agarchI(g05hkc), see above. */
  nag_generate_agarchI (num, ip, iq, &PARAM(1), gamma, &HT(1), &YT(1),
                        fcall, &RVEC(1), &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_generate_agarchI (g05hkc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  for (i = 1; i <= num; ++i)
    {
      xterm = 0.0;
      for (k = 1; k <= nreg; ++k)
        xterm += X (i, k) * BX (k);

      if (mn == 1)
        YT (i) = mean + xterm + YT (i);
      else
        YT (i) = xterm + YT (i);
    }

  for (i = 1; i <= npar; ++i)
    THETA (i) = PARAM (i) * 0.5;

  if (isym == 1)
    THETA (npar + isym) = gamma * 0.5;

  if (mn == 1)
    THETA (npar + isym + 1) = mean * 0.5;

  for (i = 1; i <= nreg; ++i)
    THETA (npar + isym + mn + i) = BX (i) * 0.5;

  maxit = 50;
  tol = 1e-12;
  stat_opt = Nag_Garch_Stationary_True;
  est_opt = Nag_Garch_Est_Initial_True;

  /* nag_estimate_agarchI(g13fac).
   * Univariate time series, parameter estimation for either a
   * symmetric GARCH process or a GARCH process with asymmetry
   * of the form (epsilon_(t~-~1)~+~gamma~)^2
   */
  nag_estimate_agarchI (&YT (1), &X (1, 1), tdx, num, ip, iq, nreg, mn, isym,
                        &THETA (1), &SE (1), &SC (1), &COVAR (1, 1), tdc, &hp,
                        &ET (1), &HT (1), &lgf, stat_opt, est_opt, maxit, tol,
                        &fail);

  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_estimate_agarchI (g13fac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf ("       Parameter estimates     Standard errors       "
           "Correct values\n");
  for (j = 1; j <= npar; ++j)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (j), SE (j),
             PARAM(j));

  if (isym)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (npar+isym),
             SE (npar+isym), gamma);
  if (mn)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (npar+isym+1),
             SE (npar+isym+1), mean);
  for (j = 1; j <= nreg; ++j)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (npar+isym+mn+j),
             SE(npar+isym+mn+j), BX(j));

  /* now forecast nt steps ahead */

  if (isym)
    {
      gamma = THETA(npar+isym);
    }
  else
    {
      gamma = 0.0;
    }

  /* nag_forecast_agarchI(g13fbc).
   * Univariate time series, forecast function for either a
   * symmetric GARCH process or a GARCH process with asymmetry
   * of the form (epsilon_(t~-~1)~+~gamma~)^2
   */
  nag_forecast_agarchI(num, nt, ip, iq, &THETA(1), gamma, &CVAR(1), &HT(1),
                       &ET(1), &fail);

  Vprintf ("\n%ld step forecast = %8.4f\n",nt,CVAR(nt));



 END:
  if (bx) NAG_FREE (bx);
  if (covar) NAG_FREE (covar);
  if (et) NAG_FREE (et);
  if (ht) NAG_FREE (ht);
  if (param) NAG_FREE (param);
  if (rvec) NAG_FREE (rvec);
  if (sc) NAG_FREE (sc);
  if (se) NAG_FREE (se);
  if (theta) NAG_FREE (theta);
  if (x) NAG_FREE (x);
  if (cvar) NAG_FREE (cvar);
  if (yt) NAG_FREE (yt);

  return exit_status;
}