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

Fortranによるサンプルソースコード : 使用ルーチン名:g05pjf

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > VARMAモデルの多変量時系列の実現値の生成

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

概要

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

VARMAモデルのデータ 

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

入力データ

(本ルーチンの詳細はg05pjf のマニュアルページを参照)
1
2
3
4
5
6
7
8
9

このデータをダウンロード
G05PJF Example Program Data
1  1  1762543     :: GENID,SUBID,SEED(1)
48 2              :: N,NREAL
2  1  0           :: K, IP, IQ
0.80  0.07        
0.00  0.58        :: End of PHI
5.00  9.00        :: XMEAN
2.97                        
0.64  5.38        :: End of VAR (lower triangle)

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に使用する生成器(genid=1:nAG基本生成器)、生成器に関する情報(subid=1:GENIDが1の場合この値は参照されません)、生成器の初期値(seed=1762543)を指定しています。
  • 3行目に生成される観測値の数(n=48)と実現値の数(nreal=2)を指定しています。
  • 3行目に多変量時系列の次数(k=2)、自己回帰パラメータ行列の数(ip=1)、移動平均パラメータ行列の数(iq=0)を指定しています。
  • 5~6行目にモデルの自己回帰パラメータ行列(phi)を指定しています。
  • 7行目に多変量時系列の平均ベクトル(xmean)を指定しています。
  • 8~9行目に分散共分散行列の下三角部分の要素(var)を指定しています。

出力結果

(本ルーチンの詳細はg05pjf のマニュアルページを参照)
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

この出力例をダウンロード
 G05PJF 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が出力されています。

ソースコード

(本ルーチンの詳細はg05pjf のマニュアルページを参照)

※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法

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

このソースコードをダウンロード
    PROGRAM g05pjfe

!      G05PJF Example Program Text

!      Mark 23 Release. nAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g05kff, g05pjf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: lseed = 1, nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: genid, i, ifail, ii, ip, iq, j, k,   &
                                          k2, l, ldvar, ldx, lphi, lr, lstate, &
                                          ltheta, mode, n, nreal, rn, subid
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: phi(:), r(:), theta(:), var(:,:),    &
                                          x(:,:), xmean(:)
       INTEGER                         :: seed(lseed)
       INTEGER, ALLOCATABLE            :: state(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          max
!      .. Executable Statements ..
       WRITE (nout,*) 'G05PJF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read in the base generator information and seed
       READ (nin,*) genid, subid, seed(1)

!      Initial call to initialiser to get size of STATE array
       lstate = 0
       ALLOCATE (state(lstate))
       ifail = 0
       CALL g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!      Reallocate STATE
       DEALLOCATE (state)
       ALLOCATE (state(lstate))

!      Initialize the generator to a repeatable sequence
       ifail = 0
       CALL g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!      Read in the sample size and number of realizations
       READ (nin,*) n, nreal

!      Read in the number of coefficients
       READ (nin,*) k, ip, iq

       k2 = k**2
       rn = max(ip,iq)
       l = k*(k+1)/2
       IF (ip>0) THEN
          l = l + (ip-1)*k2
       END IF
       IF (k>=6) THEN
          lr = (5*rn**2+1)*k2 + (4*rn+3) + 4
       ELSE
          lr = ((ip+iq)**2+1)*k2 + (4*(ip+iq)+3)*k + &
             max(k*rn*(k*rn+2),k2*(ip+iq)**2+l*(l+3)+k2*(iq+1)) + 4
       END IF
       lphi = ip*k*k
       ltheta = iq*k*k
       ldvar = k
       ldx = k
       ALLOCATE (phi(lphi),theta(ltheta),var(ldvar,k),r(lr),x(ldx,n),xmean(k))

!      Read in the AR parameters
       DO l = 1, ip
          DO i = 1, k
             ii = (l-1)*k*k + i
             READ (nin,*) (phi(ii+k*(j-1)),j=1,k)
          END DO
       END DO

!      Read in the MA parameters
       DO l = 1, iq
          DO i = 1, k
             ii = (l-1)*k*k + i
             READ (nin,*) (theta(ii+k*(j-1)),j=1,k)
          END DO
       END DO

!      Read in the means
       READ (nin,*) xmean(1:k)

!      Read in the variance / covariance matrix
       READ (nin,*) (var(i,1:i),i=1,k)

!      For the first realization we need to set up the reference vector
!      as well as generate the series
       mode = 2

!      Generate NREAL realizations
D_LP:  DO rn = 1, nreal

          ifail = 0
          CALL g05pjf(mode,n,k,xmean,ip,phi,iq,theta,var,ldvar,r,lr,state,x, &
             ldx,ifail)

!         Display the results
          WRITE (nout,99999) ' Realization Number ', rn
          DO i = 1, k
             WRITE (nout,*)
             WRITE (nout,99999) '  Series number ', i
             WRITE (nout,*) '  -------------'
             WRITE (nout,*)
             WRITE (nout,99998) x(i,1:n)
          END DO
          WRITE (nout,*)

!         For subsequent realizations we use previous reference vector
          mode = 3

       END DO D_LP

99999  FORMAT (1X,A,I0)
99998  FORMAT (8(2X,F8.3))
    END PROGRAM g05pjfe


関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks