Keyword: 自己回帰和分移動平均, ARIMA, 時系列, フィルタリング
概要
本サンプルは自己回帰和分移動平均（ARIMA）モデルによる時系列のフィルタリングを行うC言語によるサンプルプログラムです。 本サンプルは長さ296の時系列データを読み込み、またARIMAモデル(4,0,2,0,0,0,0)とARIMAフィルタリングモデル(3,0,0,0,0,0,0)を読み込み、２つの後方予測とフィルタリングされた時系列を出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_tsa_arma_filter() のExampleコードです。本サンプル及び関数の詳細情報は nag_tsa_arma_filter のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_tsa_arma_filter のマニュアルページを参照）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
| このデータをダウンロード |
nag_tsa_arma_filter (g13bac) Example Program Data
296
53.8 53.6 53.5 53.5 53.4 53.1 52.7 52.4 52.2 52.0 52.0 52.4 53.0 54.0 54.9 56.0
56.8 56.8 56.4 55.7 55.0 54.3 53.2 52.3 51.6 51.2 50.8 50.5 50.0 49.2 48.4 47.9
47.6 47.5 47.5 47.6 48.1 49.0 50.0 51.1 51.8 51.9 51.7 51.2 50.0 48.3 47.0 45.8
45.6 46.0 46.9 47.8 48.2 48.3 47.9 47.2 47.2 48.1 49.4 50.6 51.5 51.6 51.2 50.5
50.1 49.8 49.6 49.4 49.3 49.2 49.3 49.7 50.3 51.3 52.8 54.4 56.0 56.9 57.5 57.3
56.6 56.0 55.4 55.4 56.4 57.2 58.0 58.4 58.4 58.1 57.7 57.0 56.0 54.7 53.2 52.1
51.6 51.0 50.5 50.4 51.0 51.8 52.4 53.0 53.4 53.6 53.7 53.8 53.8 53.8 53.3 53.0
52.9 53.4 54.6 56.4 58.0 59.4 60.2 60.0 59.4 58.4 57.6 56.9 56.4 56.0 55.7 55.3
55.0 54.4 53.7 52.8 51.6 50.6 49.4 48.8 48.5 48.7 49.2 49.8 50.4 50.7 50.9 50.7
50.5 50.4 50.2 50.4 51.2 52.3 53.2 53.9 54.1 54.0 53.6 53.2 53.0 52.8 52.3 51.9
51.6 51.6 51.4 51.2 50.7 50.0 49.4 49.3 49.7 50.6 51.8 53.0 54.0 55.3 55.9 55.9
54.6 53.5 52.4 52.1 52.3 53.0 53.8 54.6 55.4 55.9 55.9 55.2 54.4 53.7 53.6 53.6
53.2 52.5 52.0 51.4 51.0 50.9 52.4 53.5 55.6 58.0 59.5 60.0 60.4 60.5 60.2 59.7
59.0 57.6 56.4 55.2 54.5 54.1 54.1 54.4 55.5 56.2 57.0 57.3 57.4 57.0 56.4 55.9
55.5 55.3 55.2 55.4 56.0 56.5 57.1 57.3 56.8 55.6 55.0 54.1 54.3 55.3 56.4 57.2
57.8 58.3 58.6 58.8 58.8 58.6 58.0 57.4 57.0 56.4 56.3 56.4 56.4 56.0 55.2 54.0
53.0 52.0 51.6 51.6 51.1 50.4 50.0 50.0 52.0 54.0 55.1 54.5 52.8 51.4 50.8 51.2
52.0 52.8 53.8 54.5 54.9 54.9 54.8 54.4 53.7 53.3 52.8 52.6 52.6 53.0 54.3 56.0
57.0 58.0 58.6 58.5 58.3 57.8 57.3 57.0
4 0 2 0 0 0 0
0.000
2.420 -2.380 1.160 -0.230 0.310 -0.470
3 0 0 0 0 0 0
1.970 -1.370 0.340
- １行目はタイトル行で読み飛ばされます。
- ２行目に時系列データ点の合計数(nx)を指定しています。
- ３～２１行目に時系列データ(x)を指定しています。
- ２２行目にARIMAモデルの次数ベクトル(mrx)を指定しています。
- ２３行目にARIMAモデルの定数項(cx)を指定しています。
- ２４行目にARIMAモデルの引数(parx)を指定しています。
- ２５行目にフィルタリングARIMAモデルの次数(mr)を指定しています。
- ２６行目にフィルタリングモデルの引数(par)を指定しています。
出力結果
（本関数の詳細はnag_tsa_arma_filter のマニュアルページを参照）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
| この出力例をダウンロード |
nag_tsa_arma_filter (g13bac) Example Program Results
Original Filtered
Backforecasts y-series series
-2 49.9807 3.4222
-1 52.6714 3.0809
Filtered Filtered Filtered Filtered
series series series series
1 2.9813 2 2.7803 3 3.7057 4 3.2450
5 3.0760 6 3.0070 7 3.0610 8 3.1720
9 3.1170 10 3.0360 11 3.2580 12 3.4520
13 3.3320 14 3.6980 15 3.3140 16 3.8070
17 3.3330 18 2.9580 19 3.2800 20 3.0960
21 3.2270 22 3.0830 23 2.6410 24 3.1870
25 2.9910 26 3.1110 27 2.8460 28 3.0240
29 2.7030 30 2.6130 31 2.8060 32 2.9560
33 2.8170 34 2.8950 35 2.8510 36 2.9160
37 3.2530 38 3.3050 39 3.1830 40 3.3760
41 2.9730 42 2.8610 43 3.0490 44 2.8420
45 2.3190 46 2.3660 47 2.9410 48 2.3810
49 3.3420 50 2.9340 51 3.1800 52 2.9230
53 2.6470 54 2.8860 55 2.5310 56 2.6200
57 3.4170 58 3.4940 59 3.2590 60 3.1310
61 3.1420 62 2.6710 63 2.8990 64 2.8180
65 3.2150 66 2.8800 67 2.9610 68 2.8800
69 3.0020 70 2.8930 71 3.1210 72 3.2210
73 3.2040 74 3.5360 75 3.7520 76 3.5630
77 3.7260 78 3.1560 79 3.6310 80 2.9380
81 3.1480 82 3.4490 83 3.1400 84 3.7380
85 4.1200 86 3.1540 87 3.7480 88 3.3280
89 3.3640 90 3.3400 91 3.3950 92 3.0720
93 3.0050 94 2.8520 95 2.7810 96 3.1950
97 3.2490 98 2.6370 99 3.0080 100 3.2410
101 3.5570 102 3.2080 103 3.0880 104 3.3980
105 3.1660 106 3.1960 107 3.2460 108 3.2870
109 3.1590 110 3.2620 111 2.7280 112 3.4130
113 3.2190 114 3.6750 115 3.8550 116 4.0100
117 3.5380 118 3.8440 119 3.4660 120 3.0640
121 3.4780 122 3.1140 123 3.5300 124 3.2400
125 3.3630 126 3.2610 127 3.3020 128 3.1150
129 3.3280 130 2.8730 131 3.0800 132 2.8390
133 2.6570 134 3.0260 135 2.4580 136 3.2600
137 2.8380 138 3.2150 139 3.1140 140 3.1050
141 3.1400 142 2.9100 143 3.1370 144 2.7500
145 3.1160 146 3.0680 147 2.8590 148 3.3840
149 3.5500 150 3.4160 151 3.1770 152 3.3390
153 3.0190 154 3.1780 155 3.0110 156 3.1940
157 3.2680 158 3.0500 159 2.8060 160 3.1850
161 3.0560 162 3.2690 163 2.7940 164 3.0900
165 2.7100 166 2.7890 167 2.9510 168 3.2440
169 3.2570 170 3.4360 171 3.4450 172 3.3780
173 3.3520 174 3.9180 175 2.9190 176 3.1780
177 2.2580 178 3.5150 179 2.8010 180 3.6030
181 3.2610 182 3.5300 183 3.3270 184 3.4420
185 3.5240 186 3.2720 187 3.1110 188 2.8240
189 3.2330 190 3.1500 191 3.5710 192 3.0810
193 2.7820 194 2.9040 195 3.2350 196 2.7970
197 3.1320 198 3.1680 199 4.5210 200 2.6650
201 4.6870 202 3.9470 203 3.2220 204 3.3410
205 3.9950 206 3.4820 207 3.3630 208 3.4550
209 3.2950 210 2.6910 211 3.4600 212 2.9440
213 3.4400 214 3.1830 215 3.4200 216 3.4100
217 4.0550 218 2.9990 219 3.8250 220 3.1340
221 3.5010 222 3.0430 223 3.2660 224 3.3660
225 3.2650 226 3.3720 227 3.2880 228 3.5470
229 3.6840 230 3.3100 231 3.6790 232 3.1780
233 2.9360 234 2.7910 235 3.8020 236 2.6100
237 4.1690 238 3.7460 239 3.4560 240 3.3910
241 3.5820 242 3.6220 243 3.4870 244 3.5770
245 3.4240 246 3.3960 247 3.1220 248 3.4300
249 3.4580 250 3.0280 251 3.7660 252 3.3770
253 3.2470 254 3.0180 255 2.9720 256 2.8000
257 3.2040 258 2.8020 259 3.4100 260 3.1680
261 2.4600 262 2.8810 263 3.1750 264 3.1740
265 4.8640 266 3.0600 267 2.9600 268 2.2530
269 2.5620 270 3.3150 271 3.3480 272 3.5900
273 3.2560 274 3.2320 275 3.6160 276 3.1700
277 3.2890 278 3.1200 279 3.3300 280 2.9910
281 2.9420 282 3.4070 283 2.8720 284 3.3470
285 3.1920 286 3.4880 287 4.0680 288 3.7550
289 3.0510 290 3.9680 291 3.3900 292 3.1380
293 3.6170 294 3.1700 295 3.4150 296 3.4830
- ５～６行目に元の時系列について後方予測した時系列の値とフィルタリングされた時系列の値が出力されています。
- １０～８３行目にフィルタリングされた時系列が出力されています。
ソースコード
（本関数の詳細はnag_tsa_arma_filter のマニュアルページを参照）
※本サンプルソースコードは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 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
| このソースコードをダウンロード |
/* nag_tsa_arma_filter (g13bac) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg13.h>
int main(void)
{
/* Scalars */
double a1, a2, cx, cy;
Integer i, idd, ii, ij, iqxd,
j, k, n, nb, ni, npar, nparx, nx, ny,
nser, npara, tdxxy, tdmrx, ldparx, tdparx;
Integer exit_status = 0;
/* Arrays */
double *b = 0, *fsd = 0, *fva = 0, *par = 0, *parx = 0,
*x = 0, *y = 0, *rms = 0, *parxx = 0;
Integer mr[14], mrx[7], *mrxx = 0;
Nag_ArimaOrder arimaj, arimaf, arimav;
Nag_TransfOrder transfj;
Nag_G13_Opt options;
NagError fail;
INIT_FAIL(fail);
exit_status = 0;
/* Initialize the options structure used by nag_tsa_multi_inp_model_forecast
* (g13bjc) */
/* nag_tsa_options_init (g13bxc).
* Initialization function for option setting
*/
nag_tsa_options_init(&options);
printf("nag_tsa_arma_filter (g13bac) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%ld%*[^\n] ", &nx);
if (nx > 0) {
/* Allocate array x */
if (!(x = nAG_ALLOC(nx + 2, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= nx; ++i)
scanf("%lf", &x[i - 1]);
scanf("%*[^\n] ");
/* Read univariate Arima for series */
for (i = 1; i <= 7; ++i)
scanf("%ld", &mrx[i - 1]);
scanf("%*[^\n] ");
scanf("%lf%*[^\n] ", &cx);
nparx = mrx[0] + mrx[2] + mrx[3] + mrx[5];
arimaj.p = mrx[0];
arimaj.d = mrx[1];
arimaj.q = mrx[2];
arimaj.bigp = mrx[3];
arimaj.bigd = mrx[4];
arimaj.bigq = mrx[5];
arimaj.s = mrx[6];
nser = 1;
if (nparx > 0) {
/* Allocate array parx */
if (!(parx = nAG_ALLOC(nparx + nser, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= nparx; ++i)
scanf("%lf", &parx[i - 1]);
scanf("%*[^\n] ");
/* Read model by which to filter series */
for (i = 1; i <= 7; ++i)
scanf("%ld", &mr[i - 1]);
scanf("%*[^\n] ");
arimaf.p = mr[0];
arimaf.d = mr[1];
arimaf.q = mr[2];
arimaf.bigp = mr[3];
arimaf.bigd = mr[4];
arimaf.bigq = mr[5];
arimaf.s = mr[6];
npar = mr[0] + mr[2] + mr[3] + mr[5];
if (npar > 0) {
/* Allocate array par */
if (!(par = nAG_ALLOC(npar + nparx, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= npar; ++i)
scanf("%lf", &par[i - 1]);
scanf("%*[^\n] ");
/* Initially backforecast QY values */
/* (1) Reverse series in situ */
n = nx / 2;
ni = nx;
for (i = 1; i <= n; ++i) {
a1 = x[i - 1];
a2 = x[ni - 1];
x[i - 1] = a2;
x[ni - 1] = a1;
--ni;
}
idd = mrx[1] + mrx[4];
/* (2) Possible sign reversal for ARIMA constant */
if (idd % 2 != 0)
cx = -cx;
/* (3) Calculate number of backforecasts required */
iqxd = mrx[2] + mrx[5] * mrx[6];
/* Calculate series length */
ny = nx + iqxd;
/* Allocate array y */
if (!(y = nAG_ALLOC(ny, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
if (iqxd != 0) {
/* Allocate arrays fsd, fva and st. */
if (!(fsd = nAG_ALLOC(iqxd, double)) ||
!(fva = nAG_ALLOC(iqxd, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* (4) Set up parameter list for call to forecast
* routine g13bjc
*/
npara = nparx + nser;
parx[npara - 1] = cx;
tdxxy = nser;
tdmrx = nser - 1;
ldparx = nser - 1;
tdparx = nser - 1;
if (!(rms = nAG_ALLOC(nser, double)) ||
!(parxx = nAG_ALLOC(nser, double)) ||
!(mrxx = nAG_ALLOC(7 * nser, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* nag_tsa_transf_orders (g13byc).
* Allocates memory to transfer function model orders
*/
nag_tsa_transf_orders(nser, &transfj, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_tsa_transf_orders (g13byc)"
".\n%s\n", fail.message);
exit_status = 1;
goto END;
}
rms[0] = 0;
transfj.nag_b = 0;
transfj.nag_q = 0;
transfj.nag_p = 0;
transfj.nag_r = 1;
for (i = 1; i <= 7; ++i)
mrxx[i - 1] = 0;
parxx[0] = 0;
/* Tell nag_tsa_multi_inp_model_forecast (g13bjc) not to
* print parameters on entry */
options.list = Nag_FALSE;
/* nag_tsa_multi_inp_model_forecast (g13bjc).
* Forecasting function
*/
nag_tsa_multi_inp_model_forecast(&arimaj, nser, &transfj,
parx, npara, nx, iqxd, x,
tdxxy, rms, mrxx, tdmrx,
parxx, ldparx, tdparx,
fva, fsd, &options, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_tsa_multi_inp_model_forecast "
"(g13bjc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
j = iqxd;
for (i = 1; i <= iqxd; ++i) {
y[i - 1] = fva[j - 1];
--j;
}
/* Move series into y */
j = iqxd + 1;
k = nx;
for (i = 1; i <= nx; ++i) {
if (j > 305)
goto END;
y[j - 1] = x[k - 1];
++j;
--k;
}
}
/* Move ARIMA for series into mr */
for (i = 1; i <= 7; ++i)
mr[i + 6] = mrx[i - 1];
arimav.p = mr[7];
arimav.d = mr[8];
arimav.q = mr[9];
arimav.bigp = mr[10];
arimav.bigd = mr[11];
arimav.bigq = mr[12];
arimav.s = mr[13];
/* Move parameters of ARIMA for y into par */
for (i = 1; i <= nparx; ++i)
par[npar + i - 1] = parx[i - 1];
npar += nparx;
/* Move constant and reset sign reversal */
cy = cx;
if (idd % 2 != 0)
cy = -cy;
/* Set parameters for call to filter routine
* nag_tsa_arma_filter (g13bac) */
nb = ny + MAX(mr[2] + mr[5] * mr[6],
mr[0] + mr[1] + (mr[3] + mr[4]) * mr[6]);
/* Allocate array b */
if (!(b = nAG_ALLOC(nb, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Filter series by call to nag_tsa_arma_filter (g13bac) */
/* nag_tsa_arma_filter (g13bac).
* Multivariate time series, filtering (pre-whitening) by an
* ARIMA model
*/
nag_tsa_arma_filter(y, ny, &arimaf, &arimav, par, npar, cy, b,
nb, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_tsa_arma_filter (g13bac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf(" Original Filtered\n");
printf("Backforecasts y-series series\n");
if (iqxd != 0) {
ij = -iqxd;
for (i = 1; i <= iqxd; ++i) {
printf("%8ld%17.4f%15.4f\n", ij, y[i - 1], b[i - 1]);
++ij;
}
}
printf("\n");
printf(" Filtered Filtered "
"Filtered Filtered\n");
printf(" series series "
" series series\n");
for (i = iqxd + 1; i <= ny; i += 4) {
for (ii = i; ii <= MIN(ny, i + 3); ++ii) {
printf("%5ld", ii - iqxd);
printf("%9.4f ", b[ii - 1]);
}
printf("\n");
}
}
}
}
END:
/* Free the options structure used by nag_tsa_multi_inp_model_forecast
* (g13bjc) */
/* nag_tsa_free (g13xzc).
* Freeing function for use with g13 option setting
*/
nag_tsa_free(&options);
nAG_FREE(b);
nAG_FREE(fsd);
nAG_FREE(fva);
nAG_FREE(par);
nAG_FREE(parx);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(rms);
nAG_FREE(parxx);
nAG_FREE(mrxx);
return exit_status;
}