Rでデータサイエンス:Dive into R function - vars:::Phi.varest

Rでデータサイエンス

Dive into R function - vars:::Phi.varest

Description of function

定常VARをVMA表現(vector moving average representation)するための係数行列を求める関数です。

具体的には VMA表現\[y_t=\Phi_0u_t+\Phi_1u_{t-1}+\Phi_2u_{t-2}+\cdots\]の係数行列\(\Phi\)を求めます。ここで\(u\)はイノベーション。

\(k\)を内生変数の数として\(\Phi_0\)は\(k\)次の単位行列\(I\)とすると、係数行列\(\Phi_s\)は以下の通り再帰的に求められます。

\[\Phi_s=\displaystyle\sum_{j=1}^s\Phi_{s-j}A_j\quad s=1,2,\dots\] \(A_j\)は定常VARの係数行列であり、VARモデルのラグ次数を\(p\)とすると、\(j>p\,\)の場合\(A_j\)はゼロと設定されます。

なお、\(s\)は引数 nstep の値に1を足した数値に設定されます。

Load packages

library(vars)
library(dplyr)

Package version

packageVersion("vars")

[1] '1.6.1'

Sample data

VARモデルのサンプルデータを作成します。

set.seed(20220422)
n <- 54
x0 <- rnorm(n = n)
y0 <- rep(mean(x0), 2)
for (i in 3:n) {
    y0[i] <- y0[i - 1] - y0[i - 2] + x0[i - 1] - x0[i - 2] + rnorm(1, mean = 0, sd = 0.1)
}
x <- tail(x0, -4)
y <- tail(y0, -4)
df <- data.frame(x, y)
df %>%
    head()
df %>%
    tail()

           x          y
1 -0.4561261 -2.8266550
2 -0.2478359 -0.7520638
3 -1.5430870  2.1470816
4 -0.0331243  1.6590939
5  2.0206307  1.1129648
6  1.1459143  1.5251849
            x         y
45 -1.3539026 -2.394264
46 -0.7769842  4.473228
47  0.5306799  7.431717
48 -0.4417267  4.176271
49 -0.4959663 -4.325015
50  1.1507611 -8.619026

単位根を確認します。

tseries::adf.test(x)

    Augmented Dickey-Fuller Test

data:  x
Dickey-Fuller = -4.727, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

tseries::adf.test(y)

    Augmented Dickey-Fuller Test

data:  y
Dickey-Fuller = -7.8363, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

VAR model

関数Phi.varestに渡すVARモデルを作成します。

p <- VARselect(y = df)$selection["AIC(n)"]
p

AIC(n) 
     2 

varmodel <- VAR(y = df, p = p, type = "none", ic = "AIC")
varmodel %>%
    summary(equations = "y")

VAR Estimation Results:
========================= 
Endogenous variables: x, y 
Deterministic variables: none 
Sample size: 48 
Log Likelihood: -7.759 
Roots of the characteristic polynomial:
0.9978 0.9978 0.4419 0.4419
Call:
VAR(y = df, p = p, type = "none", ic = "AIC")


Estimation results for equation y: 
================================== 
y = x.l1 + y.l1 + x.l2 + y.l2 

      Estimate Std. Error t value            Pr(>|t|)    
x.l1  0.979722   0.014740   66.47 <0.0000000000000002 ***
y.l1  0.997018   0.003474  287.00 <0.0000000000000002 ***
x.l2 -0.962741   0.015001  -64.18 <0.0000000000000002 ***
y.l2 -0.996035   0.003427 -290.62 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.0856 on 44 degrees of freedom
Multiple R-Squared: 0.9996, Adjusted R-squared: 0.9996 
F-statistic: 3.021e+04 on 4 and 44 DF,  p-value: < 0.00000000000000022 



Covariance matrix of residuals:
       x        y
x 0.7729 0.006300
y 0.0063 0.007298

Correlation matrix of residuals:
        x       y
x 1.00000 0.08388
y 0.08388 1.00000

Function code

関数Phi.varestのコードの出力をパート毎に確認します。

Phi

function(x, nstep = 10, ...) {
    UseMethod("Phi", x)
}

methods("Phi")

[1] Phi.svarest* Phi.svecest* Phi.varest*  Phi.vec2var*
see '?methods' for accessing help and source code

getS3method(f = "Phi", class = "varest")

function(x, nstep = 10, ...) {
    if (!(class(x) == "varest")) {
        stop("\nPlease provide an object of class 'varest', generated by 'VAR()'.\n")
    }
    nstep <- abs(as.integer(nstep))
    K <- x$K
    p <- x$p
    A <- as.array(Acoef(x))
    if (nstep >= p) {
        As <- array(0, dim = c(K, K, nstep + 1))
        for (i in (p + 1):(nstep + 1)) {
            As[, , i] <- matrix(0, nrow = K, ncol = K)
        }
    } else {
        As <- array(0, dim = c(K, K, p))
    }
    for (i in 1:p) {
        As[, , i] <- A[[i]]
    }
    Phi <- array(0, dim = c(K, K, nstep + 1))
    Phi[, , 1] <- diag(K)
    Phi[, , 2] <- Phi[, , 1] %*% As[, , 1]
    if (nstep > 1) {
        for (i in 3:(nstep + 1)) {
            tmp1 <- Phi[, , 1] %*% As[, , i - 1]
            tmp2 <- matrix(0, nrow = K, ncol = K)
            idx <- (i - 2):1
            for (j in 1:(i - 2)) {
                tmp2 <- tmp2 + Phi[, , j + 1] %*% As[, , idx[j]]
            }
            Phi[, , i] <- tmp1 + tmp2
        }
    }
    return(Phi)
}

x <- varmodel
nstep <- 3

6-8行目

K <- x$K
K

[1] 2

p <- x$p
p

AIC(n) 
     2 

A <- as.array(Acoef(x))
A

[[1]]
         x.l1       y.l1
x -0.02380917 0.01292143
y  0.97972182 0.99701842

[[2]]
        x.l2         y.l2
x -0.1913533  0.003956531
y -0.9627414 -0.996035353

(参考)

Bcoef(x)

         x.l1       y.l1       x.l2         y.l2
x -0.02380917 0.01292143 -0.1913533  0.003956531
y  0.97972182 0.99701842 -0.9627414 -0.996035353

10-13行目

As <- array(0, dim = c(K, K, nstep + 1))
for (i in (p + 1):(nstep + 1)) {
    As[, , i] <- matrix(0, nrow = K, ncol = K)
}
As

, , 1

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 2

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 3

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 4

     [,1] [,2]
[1,]    0    0
[2,]    0    0

17-19行目

for (i in 1:p) {
    As[, , i] <- A[[i]]
}
As

, , 1

            [,1]       [,2]
[1,] -0.02380917 0.01292143
[2,]  0.97972182 0.99701842

, , 2

           [,1]         [,2]
[1,] -0.1913533  0.003956531
[2,] -0.9627414 -0.996035353

, , 3

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 4

     [,1] [,2]
[1,]    0    0
[2,]    0    0

20-21行目

Phi <- array(0, dim = c(K, K, nstep + 1))
Phi[, , 1] <- diag(K)
Phi[, , 1]

     [,1] [,2]
[1,]    1    0
[2,]    0    1

22行目

Phi[, , 2] <- Phi[, , 1] %*% As[, , 1]
Phi[, , 2]

            [,1]       [,2]
[1,] -0.02380917 0.01292143
[2,]  0.97972182 0.99701842

24-32行目

for (i in 3:(nstep + 1)) {
    tmp1 <- Phi[, , 1] %*% As[, , i - 1]
    tmp2 <- matrix(0, nrow = K, ncol = K)
    idx <- (i - 2):1
    for (j in 1:(i - 2)) {
        tmp2 <- tmp2 + Phi[, , j + 1] %*% As[, , idx[j]]
    }
    Phi[, , i] <- tmp1 + tmp2
}
Phi

, , 1

     [,1] [,2]
[1,]    1    0
[2,]    0    1

, , 2

            [,1]       [,2]
[1,] -0.02380917 0.01292143
[2,]  0.97972182 0.99701842

, , 3

             [,1]       [,2]
[1,] -0.178127056 0.01653178
[2,] -0.009267058 0.01066978

, , 4

            [,1]         [,2]
[1,]  0.01255358  0.001216437
[2,] -1.13666989 -0.978671071

Reference

• https://cran.r-project.org/package=vars
• https://cran.r-project.org/web/packages/vars/vars.pdf#page=20
• https://www.pfaffikus.de/
• 村尾博,『Rで学ぶVAR実証分析』,オーム社,pp.249-254

最終更新

Sys.time()

[1] "2024-03-31 12:59:25 JST"

R、Quarto、Package

R.Version()$version.string

[1] "R version 4.3.3 (2024-02-29 ucrt)"

quarto::quarto_version()

[1] '1.4.542'

packageVersion(pkg = "tidyverse")

[1] '2.0.0'

著者

• アセット･マネジメント･コンサルティング株式会社 https://www.am-consulting.co.jp
• 法人概要 https://www.contact.am-consulting.co.jp/company-brochure/
• 代表取締役 三竹道雄 https://www.contact.am-consulting.co.jp/mitake-michio/
• 公開コンテンツリスト https://www.am-consulting.co.jp/amc-all-contents.html
• Charts https://www.am-consulting.co.jp/charts.html
• 法令検索 https://www.laws.am-consulting.co.jp/
• SaECaNet https://www.saecanet.com/
• twitter https://twitter.com/amc_corporation
• youtube https://www.youtube.com/user/michiomitake
• twilog https://twilog.togetter.com/amc_corporation
• github https://github.com/am-consulting