Keyword: VARMA, 多変量, 時系列, 実現値
概要
本サンプルはVARMAモデルの多変量時系列の実現値の生成を行うC言語によるサンプルプログラムです。 本サンプルは以下の式で示されるVARMAモデルに対して１つの実現値が48個の観測値をもつ２変量時系列の２つの実現値を生成し出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_rand_varma() のExampleコードです。本サンプル及び関数の詳細情報は nag_rand_varma のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_rand_varma のマニュアルページを参照）1 2 3 4 5 6 7
| このデータをダウンロード |
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
- １行目はタイトル行で読み飛ばされます。
- ２行目に多変量時系列の次数(k)、自己回帰パラメータの数(ip)、移動平均パラメータ行列の数(iq)、生成される観測値の数(n)を指定しています。
- ３～４行目にモデルの自己回帰パラメータ行列(phi)を指定しています。
- ５行目に多変量時系列の平均ベクトル(xmean)を指定しています。
- ６～７行目に分散共分散行列の要素(var)を指定しています。
出力結果
（本関数の詳細はnag_rand_varma のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
| この出力例をダウンロード |
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
- ８～１３行目に生成された１つ目の実現値の時系列１が出力されています。
- １８～２３行目に生成された１つ目の実現値の時系列２が出力されています。
- ３０～３５行目に生成された２つ目の実現値の時系列１が出力されています。
- ４０～４５行目に生成された２つ目の実現値の時系列２が出力されています。
ソースコード
（本関数の詳細はnag_rand_varma のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
| このソースコードをダウンロード |
/* 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;
}