/* nag_estimate_garchGJR (g13fec) 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)
{
  Integer exit_status=0, i, ip, iq, j, k, maxit, mn, npar, nreg, nt, num;
  Integer 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, *etm=0, fac1, gamma, hp, *ht=0, *htm=0;
  double lgf, mean, *param=0, *rvec=0, *sc=0, *se=0, *theta=0, tol, *x=0;
  double xterm, *yt=0;

  INIT_FAIL(fail);
  num = 1000;
  mn = 1;
  mean = 4.0;
  nreg = 2;
  ip = 1;
  iq = 1;
  npar = iq + ip + 1;
  nt = 6;
  tdx = nreg;
  tdc = npar+mn+nreg+1;

#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 HTM(I) htm[(I)-1]
#define HT(I) ht[(I)-1]
#define ETM(I) etm[(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 ("nag_estimate_garchGJR(g13fec) Example Program Results \n\n");


  if (!(bx = NAG_ALLOC (nreg, double))
      || !(covar = NAG_ALLOC ((npar+mn+nreg+1) * (npar+mn+nreg+1), double))
      || !(etm = NAG_ALLOC (num, double))
      || !(ht = NAG_ALLOC (num, double))
      || !(htm = 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))
      || !(cvar = NAG_ALLOC (nt, double))
      || !(x = NAG_ALLOC (num*nreg, double))
      || !(yt = NAG_ALLOC (num, double)))
    {
      Vprintf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  seed = 11;
  gamma = 0.1;
  BX (1) = 1.5;
  BX (2) = 2.5;

  for (i = 1; i <= num; ++i)
    {
      fac1 = (double) i *0.01;
      X (i, 2) = sin (fac1) * 0.7 + 0.01;
      X (i, 1) = fac1 * 0.1 + 0.5;
    }

  PARAM (1) = 0.4;
  PARAM (2) = 0.1;
  PARAM (3) = 0.7;

  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 = 200;
  /* nag_generate_garchGJR(g05hmc).
   * Univariate time series, generate n terms of an asymmetric
   * Glosten, Jagannathan and Runkle (GJR) GARCH process
   */
  nag_generate_garchGJR (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_garchGJR (g05hmc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  fcall = Nag_Garch_Fcall_False;
  /* nag_generate_garchGJR(g05hmc), see above. */
  nag_generate_garchGJR (num, ip, iq, &PARAM (1), gamma, &HT (1), &YT (1),
                         fcall, &RVEC (1), &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_generate_garchGJR (g05hmc).\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);
    }

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

  for (i = 1; i <= npar; ++i)
    THETA (i) = PARAM (i) * 0.5;
  THETA (npar + 1) = gamma * 0.5;
  if (mn == 1)
    THETA (npar + 2) = mean * 0.5;

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

  /* nag_estimate_garchGJR(g13fec).
   * Univariate time series, parameter estimation for an
   * asymmetric Glosten, Jagannathan and Runkle (GJR) GARCH
   * process
   */
  nag_estimate_garchGJR (&YT (1), &X (1, 1), tdx, num, ip, iq, nreg, mn,
                         &THETA (1), &SE (1), &SC (1), &COVAR (1, 1), tdc, &hp,
                         &ETM (1), &HTM (1), &lgf, stat_opt, est_opt, maxit,
                         tol, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_estimate_garchGJR (g13fec).\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)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA(npar+2), SE(npar+2),
             mean);
  for (j = 1; j <= nreg; ++j)
    Vprintf ("%20.4f             (%6.4f) %20.4f\n", THETA(npar+1+mn+j),
             SE(npar+1+mn+j), BX(j));

  /* now forecast nt steps ahead */

  gamma = THETA(npar+1);

  /* nag_forecast_garchGJR(g13ffc).
   * Univariate time series, forecast function for an
   * asymmetric Glosten, Jagannathan and Runkle (GJR) GARCH
   * process
   */
  nag_forecast_garchGJR(num, nt, ip, iq, &THETA(1), gamma, &CVAR(1), &HTM(1),
                        &ETM(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 (etm) NAG_FREE (etm);
  if (ht) NAG_FREE (ht);
  if (htm) NAG_FREE (htm);
  if (param) NAG_FREE (param);
  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);
  if (rvec) NAG_FREE (rvec);

  return exit_status;

}