指数平滑化モデルの時系列の実現値の生成

nag_rand_exp_smooth (g05pmc) Example Program Data
Nag_EstimateInitialValues Nag_LinearHolt
11 5 100 0.05                                      : mode,itype,n,nf,nsim,alpha
180 135 213 181 148 204 228 225 198 200 187        : y
dependent arguments for itype=Nag_LinearHolt
0.01 1.0 1.0                                       : param[0],param[1],param[2]
dependent arguments for mode=Nag_ContinueAndUpdate
11                                                 : k

nag_rand_exp_smooth (g05pmc) Example Program Results

Initial values used:
   1        168.018  
   2          3.800  

Mean Deviation     =    2.5473e+01
Absolute Deviation =    2.1233e+01

         Observed      1-Step
 Period   Values      Forecast      Residual

   1      180.000      171.818        8.182
   2      135.000      175.782      -40.782
   3      213.000      178.848       34.152
   4      181.000      183.005       -2.005
   5      148.000      186.780      -38.780
   6      204.000      189.800       14.200
   7      228.000      193.492       34.508
   8      225.000      197.732       27.268
   9      198.000      202.172       -4.172
  10      200.000      206.256       -6.256
  11      187.000      210.256      -23.256

                                            Simulated CI         Simulated CI
 Obs.  Forecast      Estimated CI        (Gaussian Errors)    (Bootstrap Errors)
 12    213.854    163.928    263.781    161.431    258.001    173.073    248.363
 13    217.685    167.748    267.622    172.660    262.100    177.311    252.638
 14    221.516    171.556    271.475    169.259    263.107    179.344    256.921
 15    225.346    175.347    275.345    180.721    272.776    183.672    260.804
 16    229.177    179.115    279.238    184.790    263.591    186.398    264.173
    95.0% CIs were produced

/* nag_rand_exp_smooth (g05pmc) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg01.h>
#include <nagg05.h>
#include <nagg13.h>

#define BLIM(I, J) blim[J*2 + I]
#define BSIM(I, J) bsim[J*nsim + I]
#define GLIM(I, J) glim[J*2 + I]
#define GSIM(I, J) gsim[J*nsim + I]

int main(void)
{
  /* Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer en, i, ival, j, k, lstate, n, nf, nsim, p, nq;
  Integer *state = 0;
  /* nAG structures */
  NagError fail;
  Nag_TailProbability tail;
  Nag_InitialValues mode;
  Nag_ExpSmoothType itype;
  /* Double scalar and array declarations */
  double ad, alpha, dv, tmp, var, z, bvar;
  double *blim = 0, *bsim = 0, *e = 0, *fse = 0, *fv = 0;
  double *glim = 0, *gsim = 0, *init = 0, *param = 0, *r = 0;
  double *res = 0, *tsim1 = 0, *tsim2 = 0, *y = 0, *yhat = 0;
  double q[2];
  /* Character scalar and array declarations */
  char smode[40], sitype[40];
  /* Choose the base generator */
  Nag_BaseRNG genid = Nag_Basic;
  Integer subid = 0;
  /* Set the seed */
  Integer seed[] = { 1762543 };
  Integer lseed = 1;

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

  printf("nag_rand_exp_smooth (g05pmc) 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;
  }

  /* Skip headings in data file */
  scanf("%*[^\n] ");
  /* Read in the initial arguments and check array sizes */
  scanf("%39s%39s%ld%ld%ld%lf%*[^\n] ", smode,
        sitype, &n, &nf, &nsim, &alpha);

  /*
   * nag_enum_name_to_value (x04nac).
   * Converts nAG enum member name to value
   */
  mode = (Nag_InitialValues) nag_enum_name_to_value(smode);
  itype = (Nag_ExpSmoothType) nag_enum_name_to_value(sitype);

  /* Allocate arrays */
  if (!(blim = nAG_ALLOC(2 * nf, double)) ||
      !(bsim = nAG_ALLOC(nsim * nf, double)) ||
      !(e = nAG_ALLOC(1, double)) ||
      !(fse = nAG_ALLOC(nf, double)) ||
      !(fv = nAG_ALLOC(nf, double)) ||
      !(glim = nAG_ALLOC(2 * nf, double)) ||
      !(gsim = nAG_ALLOC(nsim * nf, double)) ||
      !(res = nAG_ALLOC(n, double)) ||
      !(tsim1 = nAG_ALLOC(nf, double)) ||
      !(tsim2 = nAG_ALLOC(nf, double)) ||
      !(y = nAG_ALLOC(n, double)) ||
      !(yhat = nAG_ALLOC(n, double)) || !(state = nAG_ALLOC(lstate, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize 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;
  }

  for (i = 0; i < n; i++)
    scanf("%lf ", &y[i]);
  scanf("%*[^\n] ");

  /* Read in the itype dependent arguments (skipping headings) */
  scanf("%*[^\n] ");
  if (itype == Nag_SingleExponential) {
    /* Single exponential smoothing required */
    if (!(param = nAG_ALLOC(1, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf("%lf%*[^\n] ", &param[0]);
    p = 0;
    ival = 1;
  }
  else if (itype == Nag_BrownsExponential) {
    /* Browns exponential smoothing required */
    if (!(param = nAG_ALLOC(2, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf("%lf %lf%*[^\n] ", &param[0], &param[1]);
    p = 0;
    ival = 2;
  }
  else if (itype == Nag_LinearHolt) {
    /* Linear Holt smoothing required */
    if (!(param = nAG_ALLOC(3, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf("%lf %lf %lf%*[^\n] ", &param[0], &param[1], &param[2]);
    p = 0;
    ival = 2;
  }
  else if (itype == Nag_AdditiveHoltWinters) {
    /* Additive Holt Winters smoothing required */
    if (!(param = nAG_ALLOC(4, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf("%lf %lf %lf %lf %ld%*[^\n] ", &param[0], &param[1],
          &param[2], &param[3], &p);
    ival = p + 2;
  }
  else if (itype == Nag_MultiplicativeHoltWinters) {
    /* Multiplicative Holt Winters smoothing required */
    if (!(param = nAG_ALLOC(4, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf("%lf %lf %lf %lf %ld%*[^\n] ", &param[0], &param[1],
          &param[2], &param[3], &p);
    ival = p + 2;
  }
  else {
    printf("%s is an unknown type\n", sitype);
    exit_status = -1;
    goto END;
  }

  /* Allocate arrays */
  if (!(init = nAG_ALLOC(p + 2, double)) || !(r = nAG_ALLOC(p + 13, double)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in the mode dependent arguments (skipping headings) */
  scanf("%*[^\n] ");
  if (mode == Nag_InitialValuesSupplied) {
    /* User supplied initial values */
    for (i = 0; i < ival; i++)
      scanf("%lf ", &init[i]);
    scanf("%*[^\n] ");
  }
  else if (mode == Nag_ContinueAndUpdate) {
    /* Continuing from a previously saved R */
    for (i = 0; i < p + 13; i++)
      scanf("%lf ", &r[i]);
    scanf("%*[^\n] ");
  }
  else if (mode == Nag_EstimateInitialValues) {
    /* Initial values calculated from first k observations */
    scanf("%ld%*[^\n] ", &k);
  }
  else {
    printf("%s is an unknown mode\n", smode);
    exit_status = -1;
    goto END;
  }

  /* Fit a smoothing model (parameter r in
   * nag_rand_exp_smooth (g05pmc) and state in g13amc are in
   the same format) */
  nag_tsa_exp_smooth(mode, itype, p, param, n, y, k, init, nf, fv, fse, yhat,
                     res, &dv, &ad, r, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_tsa_exp_smooth (g13amc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Simulate forecast values from the model, and don't update r */
  var = dv * dv;
  en = n;

  /* Change the mode used to continue from fit model */
  mode = Nag_ContinueAndUpdate;

  /* Simulate nsim forecasts */
  for (i = 0; i < nsim; i++) {
    /* Simulations assuming Gaussian errors */
    nag_rand_exp_smooth(mode, nf, itype, p, param, init, var, r, state,
                        e, 0, tsim1, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_rand_exp_smooth (g05pmc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Bootstrapping errors */
    bvar = 0.0e0;
    nag_rand_exp_smooth(mode, nf, itype, p, param, init, bvar, r, state,
                        res, en, tsim2, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_rand_exp_smooth (g05pmc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Copy and transpose the simulated values */
    for (j = 0; j < nf; j++) {
      GSIM(i, j) = tsim1[j];
      BSIM(i, j) = tsim2[j];
    }
  }

  /* Calculate CI based on the quantiles for each simulated forecast */
  q[0] = alpha / 2.0e0;
  q[1] = 1.0e0 - q[0];
  nq = 2;
  for (i = 0; i < nf; i++) {
    nag_double_quantiles(nsim, &GSIM(0, i), nq, q, &GLIM(0, i), &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_double_quantiles (g01amc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    nag_double_quantiles(nsim, &BSIM(0, i), nq, q, &BLIM(0, i), &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_double_quantiles (g01amc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

  /* Display the forecast values and associated prediction intervals */
  printf("Initial values used:\n");
  for (i = 0; i < ival; i++)
    printf("%4ld   %12.3f  \n", i + 1, init[i]);
  printf("\n");
  printf("Mean Deviation     = %13.4e\n", dv);
  printf("Absolute Deviation = %13.4e\n\n", ad);
  printf("         Observed      1-Step\n");
  printf(" Period   Values      Forecast      Residual\n\n");
  for (i = 0; i < n; i++)
    printf("%4ld  %11.3f  %11.3f  %11.3f\n", i + 1, y[i],
           yhat[i], res[i]);
  printf("\n");
  tail = Nag_LowerTail;
  z = nag_deviates_normal(tail, q[1], &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_deviates_normal (g01fac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  printf("                                            Simulated CI"
         "         Simulated CI\n");
  printf(" Obs.  Forecast      Estimated CI        (Gaussian Errors)"
         "    (Bootstrap Errors)\n");
  for (i = 0; i < nf; i++) {
    tmp = z * fse[i];
    printf("%3ld %10.3f %10.3f %10.3f"
           " %10.3f %10.3f %10.3f %10.3f\n", n + i + 1, fv[i], fv[i] - tmp,
           fv[i] + tmp, GLIM(0, i), GLIM(1, i), BLIM(0, i), BLIM(1, i));
  }
  printf("   %5.1f%% CIs were produced\n", 100.0e0 * (1.0e0 - alpha));

END:
  nAG_FREE(blim);
  nAG_FREE(bsim);
  nAG_FREE(e);
  nAG_FREE(fse);
  nAG_FREE(fv);
  nAG_FREE(glim);
  nAG_FREE(gsim);
  nAG_FREE(init);
  nAG_FREE(param);
  nAG_FREE(r);
  nAG_FREE(res);
  nAG_FREE(tsim1);
  nAG_FREE(tsim2);
  nAG_FREE(y);
  nAG_FREE(yhat);
  nAG_FREE(state);

  return exit_status;
}