関連情報
ホーム > 製品 > NAG数値計算ライブラリ > サンプルソースコード集 > VARMAモデルの多変量時系列の実現値の生成 (C言語/C++)

VARMAモデルの多変量時系列の実現値の生成

C言語によるサンプルソースコード
使用関数名:nag_rand_varma (g05pjc)

Keyword: VARMA, 多変量, 時系列, 実現値

概要

本サンプルはVARMAモデルの多変量時系列の実現値の生成を行うC言語によるサンプルプログラムです。 本サンプルは以下の式で示されるVARMAモデルに対して1つの実現値が48個の観測値をもつ2変量時系列の2つの実現値を生成し出力します。

VARMAモデルのデータ 

※本サンプルはNAG Cライブラリに含まれる関数 nag_rand_varma() のExampleコードです。本サンプル及び関数の詳細情報は nag_rand_varma のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はnag_rand_varma のマニュアルページを参照)

このデータをダウンロード
nag_rand_varma (g05pjc) Example Program Data
 2 1 0 48                    : k, ip, iq, n
 0.80 0.07
 0.00 0.58                   : phi(,,1)
 5.00 9.00                   : xmean
 2.97
 0.64 5.38                   : var

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に多変量時系列の次数(k)、自己回帰パラメータの数(ip)、移動平均パラメータ行列の数(iq)、生成される観測値の数(n)を指定しています。
  • 3〜4行目にモデルの自己回帰パラメータ行列(phi)を指定しています。
  • 5行目に多変量時系列の平均ベクトル(xmean)を指定しています。
  • 6〜7行目に分散共分散行列の要素(var)を指定しています。

出力結果

(本関数の詳細はnag_rand_varma のマニュアルページを参照)

この出力例をダウンロード
nag_rand_varma (g05pjc) Example Program Results

  Realization Number 1

   Series number   1
   -------------

     4.833     2.813     3.224     3.825     1.023     1.415     2.184     3.005
     5.547     4.832     4.705     5.484     9.407    10.335     8.495     7.478
     6.373     6.692     6.698     6.976     6.200     4.458     2.520     3.517
     3.054     5.439     5.699     7.136     5.750     8.497     9.563    11.604
     9.020    10.063     7.976     5.927     4.992     4.222     3.982     7.107
     3.554     7.045     7.025     4.106     5.106     5.954     8.026     7.212
 
   Series number   2
   -------------

     8.458     9.140    10.866    10.975     9.245     5.054     5.023    12.486
    10.534    10.590    11.376     8.793    14.445    13.237    11.030     8.405
     7.187     8.291     5.920     9.390    10.055     6.222     7.751    10.604
    12.441    10.664    10.960     8.022    10.073    12.870    12.665    14.064
    11.867    12.894    10.546    12.754     8.594     9.042    12.029    12.557
     9.746     5.487     5.500     8.629     9.723     8.632     6.383    12.484
 
  Realization Number 2

   Series number   1
   -------------

     5.396     4.811     2.685     5.824     2.449     3.563     5.663     6.209
     3.130     4.308     4.333     4.903     1.770     1.278     1.340    -0.527
     1.745     3.211     4.478     5.170     5.365     4.852     6.080     6.464
     2.765     2.148     6.641     7.224    10.316     7.102     5.604     3.934
     4.839     3.698     5.210     5.384     7.652     7.315     7.332     7.561
     7.537     7.788     6.868     7.575     6.108     6.188     8.132    10.310
 
   Series number   2
   -------------

    11.345    10.070    13.654    12.409    11.329    13.054    12.465     9.867
    10.263    13.394    10.553    10.331     7.814     8.747    10.025    11.167
    10.626     9.366     9.607     9.662    10.492    10.766    11.512    10.813
    10.799     8.780     9.221    14.245    11.575    10.620     8.282     5.447
     9.935     9.386    11.627    10.066    11.394     7.951     7.907    12.616
    15.246     9.962    13.216    11.350    11.227     6.021     6.968    12.428
 

  • 8〜13行目に生成された1つ目の実現値の時系列1が出力されています。
  • 18〜23行目に生成された1つ目の実現値の時系列2が出力されています。
  • 30〜35行目に生成された2つ目の実現値の時系列1が出力されています。
  • 40〜45行目に生成された2つ目の実現値の時系列2が出力されています。

ソースコード

(本関数の詳細はnag_rand_varma のマニュアルページを参照)

※本サンプルソースコードはNAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
/* nag_rand_varma (g05pjc) 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 <nagg05.h>

#define VAR(I, J)      var[(order == Nag_RowMajor)?(I*pdvar+J):(J*pdvar+I)]
#define X(I, J)        x[(order == Nag_RowMajor)?(I*pdx+J):(J*pdx+I)]
#define PHI(i, j, l)   phi[l*k*k+j*k+i]
#define THETA(i, j, l) theta[l*k*k+j*k+i]

int main(void)
{
  /* Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer lr, x_size, var_size, i, ip, iq, j, k, l, lstate, n, tmp1,
         tmp2, tmp3, tmp4, tmp5;
  Integer *state = 0;
  Integer pdx, pdvar;
  /* NAG structures */
  NagError fail;
  Nag_ModeRNG mode;
  /* Double scalar and array declarations */
  double *phi = 0, *r = 0, *theta = 0, *var = 0, *x = 0, *xmean = 0;
  /* Use column major order */
  Nag_OrderType order = Nag_ColMajor;
  /* 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_varma (g05pjc) 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;
  }

  /* Read data from a file */
  /* Skip heading */
  scanf("%*[^\n] ");
  /* Read in initial parameters */
  scanf("%ld%ld%ld%ld%*[^\n] ", &k,
        &ip, &iq, &n);

  /* Calculate the size of the reference vector */
  tmp1 = (ip > iq) ? ip : iq;
  if (ip == 0) {
    tmp2 = k * (k + 1) / 2;
  }
  else {
    tmp2 = k * (k + 1) / 2 + (ip - 1) * k * k;
  }
  tmp3 = ip + iq;
  if (k >= 6) {
    lr = (5 * tmp1 * tmp1 + 1) * k * k + (4 * tmp1 + 3) * k + 4;
  }
  else {
    tmp4 = k * tmp1 * (k * tmp1 + 2);
    tmp5 = k * k * tmp3 * tmp3 + tmp2 * (tmp2 + 3) + k * k * (iq + 1);
    lr = (tmp3 * tmp3 + 1) * k * k + (4 * tmp3 + 3) * k +
           ((tmp4 > tmp5) ? tmp4 : tmp5) + 4;
  }

  pdvar = k;
  pdx = (order == Nag_ColMajor) ? k : n;
  x_size = (order == Nag_ColMajor) ? pdx * n : pdx * k;
  var_size = pdvar * k;

  /* Allocate arrays */
  if (!(phi = NAG_ALLOC(ip * k * k, double)) ||
      !(r = NAG_ALLOC(lr, double)) ||
      !(theta = NAG_ALLOC(iq * k * k, double)) ||
      !(var = NAG_ALLOC(var_size, double)) ||
      !(x = NAG_ALLOC(x_size, double)) ||
      !(xmean = NAG_ALLOC(k, double)) ||
      !(state = NAG_ALLOC(lstate, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in the AR parameters */
  for (l = 0; l < ip; l++) {
    for (i = 0; i < k; i++) {
      for (j = 0; j < k; j++)
        scanf("%lf ", &PHI(i, j, l));
    }
  }
  scanf("%*[^\n] ");

  /* Read in the MA parameters */
  if (iq > 0) {
    for (l = 0; l < iq; l++) {
      for (i = 0; i < k; i++) {
        for (j = 0; j < k; j++)
          scanf("%lf ", &THETA(i, j, l));
      }
    }
    scanf("%*[^\n] ");
  }

  /* Read in the means */
  for (i = 0; i < k; i++)
    scanf("%lf ", &xmean[i]);
  scanf("%*[^\n] ");

  /* Read in the variance / covariance matrix */
  for (i = 0; i < k; i++) {
    for (j = 0; j <= i; j++)
      scanf("%lf ", &VAR(i, j));
  }
  scanf("%*[^\n] ");

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

  /* Generate the first realization */
  mode = Nag_InitializeAndGenerate;
  nag_rand_varma(order, mode, n, k, xmean, ip, phi, iq, theta,
                 var, pdvar, r, lr, state, x, pdx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_varma (g05pjc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Display the results */
  printf("  Realization Number 1\n");
  for (i = 0; i < k; i++) {
    printf("\n   Series number %3ld\n", i + 1);
    printf("   -------------\n\n ");
    for (j = 0; j < n; j++)
      printf("%9.3f%s", X(i, j), (j + 1) % 8 ? " " : "\n ");
  }
  printf("\n");

  /* Generate a second realization */
  mode = Nag_ReGenerateFromReference;
  nag_rand_varma(order, mode, n, k, xmean, ip, phi, iq, theta,
                 var, pdvar, r, lr, state, x, pdx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_rand_varma (g05pjc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Display the results */
  printf("  Realization Number 2\n");
  for (i = 0; i < k; i++) {
    printf("\n   Series number %3ld\n", i + 1);
    printf("   -------------\n\n ");
    for (j = 0; j < n; j++)
      printf("%9.3f%s", X(i, j), (j + 1) % 8 ? " " : "\n ");
  }
  printf("\n");

END:
  NAG_FREE(phi);
  NAG_FREE(r);
  NAG_FREE(theta);
  NAG_FREE(var);
  NAG_FREE(x);
  NAG_FREE(xmean);
  NAG_FREE(state);

  return exit_status;
}


Results matter. Trust NAG.

Privacy Policy | Trademarks