/* 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>

#define X(I, J) x[(I) *tdx + (J)]

int main(void)
{
  /* Integer scalar and array declarations */
  Integer                    exit_status = 0;
  Integer                    i, j, k, npar, tdc, tdx, lstate, lr;
  Integer                    *state = 0;

  /* NAG structures and data types */
  NagError                   fail;
  Nag_Boolean                fcall;

  /* Double scalar and array declarations */
  double                     *covar = 0, *cvar = 0, *et = 0, *ht = 0, *r = 0;
  double                     *sc = 0, *se = 0, *theta = 0, *x = 0, *yt = 0;
  double                     fac1, hp, lgf, xterm;

  /* Choose the base generator */
  Nag_BaseRNG                genid = Nag_Basic;
  Integer                    subid = 0;

  /* Set the seed */
  Integer                    seed[] = { 1762543 };
  Integer                    lseed = 1;

  /* Set parameters for the (randomly generated) time series ... */
  /* Generate data assuming normally distributed errors */
  Nag_ErrorDistn             dist = Nag_NormalDistn;
  double                     df = 0;

  /* Size of the time series */
  Integer                    num = 1000;

  /* MA and AR parameters */
  Integer                    ip = 1;
  Integer                    iq = 1;
  double                     param[] = { 0.15, 0.1, 0.8, 0.1 };

  /* Asymmetry parameter */
  double                     gamma = -0.3;

  /* Regression parameters */
  Integer                    nreg = 2;
  double                     mean = 3.0;
  double                     bx[] = { 1.5, 2.5 };
  /* ... end of parameters for (randomly generated) time series */

  /* When fitting a model to the time series ... */
  /* Include asymmetry parameter in the model */
  Integer                    isym = 1;

  /* Include mean in the model */
  Integer                    mn = 1;

  /* Use the following maaximum number of iterations and tolerance */
  Integer                    maxit = 50;
  double                     tol = 1e-12;

  /* Enforce stationary conditions */
  Nag_Garch_Stationary_Type  stat_opt = Nag_Garch_Stationary_True;

  /* Estimate initial values for regression parameters */
  Nag_Garch_Est_Initial_Type est_opt = Nag_Garch_Est_Initial_True;

  /* Set the number of values to forecast from the fitted model */
  Integer                    nt = 3;
  /* ... end of model fitting options */


  /* Initialise the error structure */
  INIT_FAIL(fail);

  printf("nag_estimate_agarchI (g13fac) Example Program Results \n\n");

  /* Get the length of the state array */
  lstate = -1;
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Derive various amounts */
  npar = iq + ip + 1;
  tdc = npar + mn + isym + nreg;
  tdx = nreg;

  /* Calculate the size of the reference vector */
  lr = 2 * (iq + ip + 2);

  /* Allocate arrays */
  if (!(covar = NAG_ALLOC((npar + mn + isym + nreg) * tdc, double)) ||
      !(et = NAG_ALLOC(num, double)) ||
      !(ht = NAG_ALLOC(num, 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)) ||
      !(state = NAG_ALLOC(lstate, Integer)) ||
      !(r = NAG_ALLOC(lr, double)) ||
      !(x = NAG_ALLOC(num * tdx, double)) ||
      !(cvar = NAG_ALLOC(nt, double)) ||
      !(yt = NAG_ALLOC(num, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Initialise the generator to a repeatable sequence */
  nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Set up the time dependent exogenous matrix x */
  for (i = 0; i < num; ++i)
    {
      fac1 = (double)(i + 1) * 0.01;
      X(i, 0) = sin(fac1) * 0.7 + 0.01;
      X(i, 1) = fac1 * 0.1 + 0.5;
    }

  /* Generate a realization of a random AGARCH I time series and discard it */
  fcall = Nag_TRUE;
  nag_rand_agarchI(dist, num, ip, iq, param, gamma, df, ht, yt, fcall, r, lr,
                   state, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_agarchI (g05pdc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Generate a realization of a random AGARCH I time series to use */
  fcall = Nag_FALSE;
  nag_rand_agarchI(dist, num, ip, iq, param, gamma, df, ht, yt, fcall, r, lr,
                   state, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_rand_agarchI (g05pdc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Adjust the randomly generated time series to take into account for the
     exogenous matrix x */
  for (i = 0; i < num; ++i)
    {
      xterm = 0.0;
      for (k = 0; k < nreg; ++k)
        xterm += X(i, k) * bx[k];

      if (mn == 1)
        yt[i] = mean + xterm + yt[i];
      else
        yt[i] = xterm + yt[i];
    }

  /* Set initial estimates for the parameters */
  for (i = 0; i < npar; ++i)
    theta[i] = param[i] * 0.5;
  if (isym == 1)
    theta[npar + isym - 1] = gamma * 0.5;
  if (mn == 1)
    theta[npar + isym] = mean * 0.5;
  for (i = 0; i < nreg; ++i)
    theta[npar + isym + mn + i] = bx[i] * 0.5;

  /* 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, x, tdx, num, ip, iq, nreg, mn, isym, theta, se, sc,
                       covar, tdc, &hp, et, ht, &lgf, stat_opt, est_opt, maxit,
                       tol, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_estimate_agarchI (g13fac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* Display the results */
  printf("       Parameter estimates     Standard errors       "
          "Correct values\n");
  for (j = 0; j < npar; ++j)
    printf("%20.4f             (%6.4f) %20.4f\n", theta[j], se[j],
            param[j]);

  if (isym)
    printf("%20.4f             (%6.4f) %20.4f\n",
            theta[npar + isym - 1], se[npar + isym - 1], gamma);
  if (mn)
    printf("%20.4f             (%6.4f) %20.4f\n", theta[npar + isym],
            se[npar + isym], mean);
  for (j = 0; j < nreg; ++j)
    printf("%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 - 1];
    }
  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, gamma, cvar, ht, et, &fail);
  printf("\n%ld step forecast = %8.4f\n", nt, cvar[nt - 1]);


 END:
  if (covar) NAG_FREE(covar);
  if (et) NAG_FREE(et);
  if (ht) NAG_FREE(ht);
  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);
  if (state) NAG_FREE(state);
  if (r) NAG_FREE(r);

  return exit_status;
}