/* nag_estimate_agarchII (g13fcc) Example Program. 
 *
 * Copyright 2000 Numerical Algorithms Group.
 *
 * NAG C Library
 *
 * Mark 6, 2000.
 *
 */

#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)
{
  double *bx=0, *covar=0, *et=0, fac1, gamma, hp;
  double *ht=0, lgf, mean, *param=0, *rvec=0;
  double *sc=0, *se=0, *theta=0, *cvar=0,tol;
  double *x=0, xterm, *yt=0;
  Integer exit_status = 0;
  Integer i, nt, ip, iq, j, k;
  Integer tdx, tdc;
  Integer maxit, mn, npar, nreg, seed, num, num_startup;
  Nag_Garch_Fcall_Type fcall;
  Nag_Garch_Stationary_Type stat_opt;
  Nag_Garch_Est_Initial_Type est_opt;
  NagError fail;
  
  INIT_FAIL(fail);
  mn = 1;
  nreg = 2;
  gamma = -0.4;
  ip = 1;
  iq = 1;
  num = 1000;
  nt = 3;
  npar = iq + ip + 1;
  
  tdc = npar+mn+nreg+1; 
  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)]
  
  
  
  Vprintf ("g13fcc Example Program Results\n\n");
  
  if (!(bx = NAG_ALLOC (nreg, double))
      || !(covar = NAG_ALLOC ((npar+mn+nreg+1) * (npar+mn+nreg+1), double))
      || !(et = NAG_ALLOC (num, double))
      || !(ht = NAG_ALLOC (num, double))
      || !(param = NAG_ALLOC (npar+mn+nreg+1, double))
      || !(rvec = NAG_ALLOC (40, double))
      || !(sc = NAG_ALLOC (npar+mn+nreg+1, double))
      || !(se = NAG_ALLOC (npar+mn+nreg+1, double))
      || !(theta = NAG_ALLOC (npar+mn+nreg+1, 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;
    }
  
  
  seed = 11;
  BX (1) = 1.5;
  BX (2) = 2.5;
  mean = 3.0;
  
  for (i = 1; i <= num; ++i)
    {
      fac1 = (double) i *.01;
      X (i, 1) = sin (fac1) * 0.7 + 0.01;
      X (i, 2) = fac1 * 0.1 + 0.5;
    }
  
  PARAM (1) = 0.2;
  PARAM (2) = 0.2;
  PARAM (3) = 0.7;
  
  fcall = Nag_Garch_Fcall_True;
  g05cbc(seed);
  num_startup = 300;
  g05hlc (num_startup, ip, iq, &PARAM (1), gamma, &HT (1), &YT (1), 
          fcall, &RVEC (1), &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from g05hlc.\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  
  fcall = Nag_Garch_Fcall_False;
  g05hlc (num, ip, iq, &PARAM (1), gamma, &HT (1), &YT (1), 
          fcall, &RVEC (1), &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from g05hlc.\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;
  
  THETA (npar + 1) = gamma * 0.5;
  if (mn == 1)
    THETA (npar + mn + 1) = mean * 0.5;
  
  for (i = 1; i <= nreg; ++i)
    THETA (npar + mn + 1 + i) = BX (i) * 0.5;
  
  maxit = 50;
  tol = 1e-12;
  
  stat_opt = Nag_Garch_Stationary_True;
  est_opt  = Nag_Garch_Est_Initial_True;
  
  g13fcc (&YT (1), &X (1, 1), tdx, num, ip, iq, nreg, mn,
          &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 g13fcc.\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));
  Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (npar+1), SE (npar+1), gamma);
  if (mn == 1)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (npar+mn+1), SE (npar+mn+1), mean);
  for (j = 1; j <= nreg; ++j) 
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA (mn+npar+1+j), SE(mn+npar+1+j), BX(j));
  

/* now forecast nt steps ahead */

  gamma = THETA(npar+1);

  g13fdc(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 (cvar) NAG_FREE (cvar);
  if (x) NAG_FREE (x);
  if (yt) NAG_FREE (yt);
  
  return exit_status;
}