Keyword: ML法, 自己回帰移動平均, VARMA, 時系列, フィット
概要
本サンプルはML法を用いた自己回帰移動平均(VARMA)モデルのフィットを行うC言語によるサンプルプログラムです。 本サンプルは以下の式で表されるARモデルを長さ48の二つの時系列にフィットさせ、結果を出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_tsa_varma_estimate() のExampleコードです。本サンプル及び関数の詳細情報は nag_tsa_varma_estimate のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
(本関数の詳細はnag_tsa_varma_estimate のマニュアルページを参照)1 2 3 4 5 6 7 8 9 10 11 12 13 14
このデータをダウンロード |
nag_tsa_varma_estimate (g13ddc) Example Program Data 2 1 0 48 Nag_MeanInclude : k, ip, iq, n, mean -1.490 -1.620 5.200 6.230 6.210 5.860 4.090 3.180 2.620 1.490 1.170 0.850 -0.350 0.240 2.440 2.580 2.040 0.400 2.260 3.340 5.090 5.000 4.780 4.110 3.450 1.650 1.290 4.090 6.320 7.500 3.890 1.580 5.210 5.250 4.930 7.380 5.870 5.810 9.680 9.070 7.290 7.840 7.550 7.320 7.970 7.760 7.000 8.350 7.340 6.350 6.960 8.540 6.620 4.970 4.550 4.810 4.750 4.760 10.880 10.010 11.620 10.360 6.400 6.240 7.930 4.040 3.730 5.600 5.350 6.810 8.270 7.680 6.650 6.080 10.250 9.140 17.750 13.300 9.630 6.800 4.080 5.060 4.940 6.650 7.940 10.760 11.890 5.850 9.010 7.500 10.020 10.380 8.150 8.370 10.730 12.140 : w
- 1行目はタイトル行で読み飛ばされます。
- 2行目に観測された時系列の数(k)、ARパラメータ行列の数(ip)、MAパラメータ行列の数(iq)、時系列の観測値の数(n)、平均μ が推定されたかゼロとして扱われるか(mean)を指定しています。 "Nag_MeanInclude"はμ が推定されたことを意味します。
- 3~14行目に時系列データWを指定しています。
出力結果
(本関数の詳細はnag_tsa_varma_estimate のマニュアルページを参照)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
この出力例をダウンロード |
nag_tsa_varma_estimate (g13ddc) Example Program Results VALUE OF LOG LIKELIHOOD FUNCTION ON EXIT = -0.20280E+03 MAXIMUM LIKELIHOOD ESTIMATES OF AR PARAMETER MATRICES ----------------------------------------------------- PHI(1) ROW-WISE : 0.802 0.065 ( 0.091)( 0.102) 0.000 0.575 ( 0.000)( 0.121) MAXIMUM LIKELIHOOD ESTIMATE OF PROCESS MEAN ------------------------------------------- 4.271 7.825 ( 1.219)( 0.776) MAXIMUM LIKELIHOOD ESTIMATE OF SIGMA MATRIX ------------------------------------------- 2.964 0.637 5.380 RESIDUAL SERIES NUMBER 1 ------------------------- T 1 2 3 4 5 6 7 8 V(T) -3.33 -1.24 5.75 1.27 0.32 0.11 -1.27 -0.73 T 9 10 11 12 13 14 15 16 V(T) -0.58 -1.26 -0.67 -1.13 -2.02 -0.57 1.24 -0.13 T 17 18 19 20 21 22 23 24 V(T) -0.77 -2.09 1.34 0.95 1.71 0.23 -0.01 -0.60 T 25 26 27 28 29 30 31 32 V(T) -0.68 -1.89 -0.77 2.05 2.11 0.94 -3.32 -2.50 T 33 34 35 36 37 38 39 40 V(T) 3.16 0.47 0.05 2.77 -0.82 0.25 3.99 0.20 T 41 42 43 44 45 46 47 48 V(T) -0.70 1.07 0.44 0.28 1.09 0.50 -0.10 1.70 RESIDUAL SERIES NUMBER 2 ------------------------- T 1 2 3 4 5 6 7 8 V(T) -0.19 -1.20 -0.02 1.21 -1.62 -2.16 -1.63 -1.13 T 9 10 11 12 13 14 15 16 V(T) -1.34 -1.30 4.82 0.43 2.54 0.35 -2.88 -0.77 T 17 18 19 20 21 22 23 24 V(T) 1.02 -3.85 -1.92 0.13 -1.20 0.41 1.03 -0.40 T 25 26 27 28 29 30 31 32 V(T) -1.09 -1.07 3.43 -0.08 9.17 -0.23 -1.34 -2.06 T 33 34 35 36 37 38 39 40 V(T) -3.16 -0.61 -1.30 0.48 0.79 2.87 2.38 -4.31 T 41 42 43 44 45 46 47 48 V(T) 2.32 -1.01 2.38 1.29 -1.14 0.36 2.59 2.64
- 4行目に終了時の対数尤度関数の値が出力されています。
- 9~13行目にARパラメータ行列の最尤推定値が出力されています。
- 18~19行目に処理平均の最尤推定値が出力されています。
- 24~26行目にシグマ行列の最尤推定値が出力されています。
- 31~47行目に一つ目の残差時系列が出力されています。
- 52~68行目に二つ目の残差時系列が出力されています。
ソースコード
(本関数の詳細はnag_tsa_varma_estimate のマニュアルページを参照)
※本サンプルソースコードは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
このソースコードをダウンロード |
/* nag_tsa_varma_estimate (g13ddc) Example Program. * * CLL6I261D/CLL6I261DL Version. * * Copyright 2017 Numerical Algorithms Group. * * Mark 26.1, 2017. */ /* Pre-processor includes */ #include <stdio.h> #include <string.h> #include <math.h> #include <nag.h> #include <nag_stdlib.h> #include <nagg13.h> int main(void) { /*Logical scalar and array declarations */ Nag_Boolean exact; Nag_IncludeMean mean; Nag_Boolean *parhld = 0; /*Integer scalar and array declarations */ Integer exit_status = 0; Integer i, ip, iprint, iq, ishow, j, k, kmax, maxcal, n; Integer niter, npar, pdcm, pdqq, pdw; /*Double scalar and array declarations */ double cgetol, rlogl; double *cm = 0, *g = 0, *par = 0, *qq = 0, *v = 0, *w = 0; /*Character scalar and array declarations */ char smean[40]; /*nAG Types */ NagError fail; INIT_FAIL(fail); printf("nag_tsa_varma_estimate (g13ddc) Example Program Results\n"); /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%ld%ld%ld%ld%39s%*[^\n] ", &k, &ip, &iq, &n, smean); npar = (ip + iq) * k * k; kmax = 3; /* * nag_enum_name_to_value (x04nac). * Converts nAG enum member name to value */ mean = (Nag_IncludeMean) nag_enum_name_to_value(smean); if (mean == Nag_MeanInclude) npar = npar + k; printf("\n"); pdcm = npar; pdqq = kmax; #define QQ(I, J) qq[(J-1)*pdqq + I-1] pdw = kmax; #define W(I, J) w[(J-1)*pdw + I-1] if (!(cm = nAG_ALLOC(npar * npar, double)) || !(g = nAG_ALLOC(npar, double)) || !(par = nAG_ALLOC(npar, double)) || !(qq = nAG_ALLOC(kmax * k, double)) || !(v = nAG_ALLOC(kmax * n, double)) || !(w = nAG_ALLOC(kmax * n, double)) || !(parhld = nAG_ALLOC(npar, Nag_Boolean))) { printf("Allocation failure\n"); exit_status = -1; goto END; } for (i = 0; i < npar; i++) { par[i] = 0.00e0; parhld[i] = Nag_FALSE; } /* Set all elements of Q to zero to use the covariance matrix */ /* between the K time series as the initial estimate of the */ /* covariance matrix */ for (j = 1; j <= k; j++) { for (i = j; i <= k; i++) QQ(i, j) = 0.00e0; } for (i = 1; i <= k; i++) { for (j = 1; j <= n; j++) scanf("%lf ", &W(i, j)); } scanf("%*[^\n] "); parhld[2] = Nag_TRUE; exact = Nag_TRUE; iprint = (-(1)); cgetol = 0.00010e0; maxcal = 40 * npar * (npar + 5); ishow = 2; /* * nag_tsa_varma_estimate (g13ddc) * Multivariate time series, estimation of varma model */ fflush(stdout); nag_tsa_varma_estimate(k, n, ip, iq, mean, par, npar, qq, kmax, w, parhld, exact, iprint, cgetol, maxcal, ishow, 0, &niter, &rlogl, v, g, cm, pdcm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_tsa_varma_estimate (g13ddc).\n%s\n", fail.message); exit_status = 1; goto END; } END: nAG_FREE(cm); nAG_FREE(g); nAG_FREE(par); nAG_FREE(qq); nAG_FREE(v); nAG_FREE(w); nAG_FREE(parhld); return exit_status; }