ブログトップ 記事一覧 ログイン 無料ブログ開設

hamadakoichi blog このページをアンテナに追加 RSSフィード Twitter

2011-04-27

[] R言語による多変量時系列分析−Dynamic Conditional Correlation GARCH Model  R言語による多変量時系列分析−Dynamic Conditional Correlation GARCH Modelを含むブックマーク  R言語による多変量時系列分析−Dynamic Conditional Correlation GARCH Modelのブックマークコメント

R言語による多変量時系列分析。

複数グループ・複数項目の時系列データで時系列間の関係性・影響を明らかにする。


以下、実行可能なR言語ソースコードを用い紹介する。


例:複数時系列間の関係性・影響
  • 農家ごとの農地の肥沃度・人の各仕事の労働量・各肥料の投入量の各時系列、および、生産量の時系列から、肥沃度・各労働量・各肥料量が生産量に与える影響を明らかにする。
  • 国ごとの各政策種類への投資額の時系列、および、GDP から、各政策投資GDPに与える影響を明らかにする
  • 地域ごとの各キャンペーン種類の投入量の各時系列、および、商品売上の時系列から、各地域の各キャンペーンが売上に与える影響を明らかにする。等。

R言語による時系列分析入門

R言語でのAR、ARMA/ARIMA、ARFIMA、ARCH、などの時系列モデルの入門資料(※1年前のTokyo.R 講師資料)

第4回R勉強会@東京 で話してきた -「R言語による時系列分析」 - hamadakoichi blog

  • AGENDA
    • 自己紹介
    • 時系列分析とは
    • データ操作
    • モデル
      • 確率密度
      • 期待値と分散
      • ARモデル
      • ARMA/ARIMAモデル
      • ARFIMAモデル
      • ARCHモデル
    • 最後に

R Library - Dynamic Conditional Correlation GARCH

複数時系列の時系列分析として、R言語では plm 、ccgarch 等のパッケージがあるが、今回は ccgarch (Conditional Correlation GARCH) を紹介する。

※GARCH: General AutoRegressive Conditional Heteroscedasticity


Dynamic Conditional Correlation GARCH (DCC- GARCH) Model 資料:

Vector Autoregressive (VAR) Model 資料:


Panel Linear Model(plm)を用いた実装は以下に記載している:

以下に、解説も含めた実装を記載する。


データ

Rの Panel Linear Model (plm) Package での組み込みデータ「Produc: 米国の48州の1970-1986のパネルデータ」を用いてみる。

library("plm")

# 米国の48州の1970-1986のパネルデータ
data("Produc", package = "plm")

# ALABAMA州の1970-1986データ
alabama <- Produc[Produc$state == "ALABAMA", ]

# 時系列オブジェクトの作成
library(zoo) # 
(d.index <- alabama$year)
data <- zoo(alabama[,3:10], d.index)

# 表示
par(cex.axis = 1.5, cex.lab = 2)
plot(data, main = "ALABAMA州データ(1970-1986)")

f:id:hamadakoichi:20110428033522j:image:w550

> head(data, 10)
         pcap     hwy   water    util       pc   gsp    emp unemp
1970 15032.67 7325.80 1655.68 6051.20 35793.80 28418 1010.5   4.7
1971 15501.94 7525.94 1721.02 6254.98 37299.91 29375 1021.9   5.2
1972 15972.41 7765.42 1764.75 6442.23 38670.30 31303 1072.3   4.7
1973 16406.26 7907.66 1742.41 6756.19 40084.01 33430 1135.5   3.9
1974 16762.67 8025.52 1734.85 7002.29 42057.31 33749 1169.8   5.5
1975 17316.26 8158.23 1752.27 7405.76 43971.71 33604 1155.4   7.7
1976 17732.86 8228.19 1799.74 7704.93 50221.57 35764 1207.0   6.8
1977 18111.93 8365.67 1845.11 7901.15 51084.99 37463 1269.2   7.4
1978 18479.74 8510.64 1960.51 8008.59 52604.05 39964 1336.5   6.3
1979 18881.49 8640.61 2081.91 8158.97 54525.86 40979 1362.0   7.1

RでのDynamic Conditional Correlation GARCH Modelを用いた実装

# 必要な R Package の Load 
library(ccgarch) # Conditional Correlation GARCH Model
library(vars) # Vector Auto Regressive Model 

(nobs <- dim(data)[1]) #サンプル数 T
(ndim <- dim(data)[2]) #サンプルの次元 N

# VARselect(y, lag.max = 10, type = c("const", "trend", "both", "none"),
#            season = NULL, exogen = NULL)
# 
# VAR(p)(ベクトル自己回帰モデル)の次数pを決定するvarパッケージの関数。
# 1次からlag.maxで指定した次数までのVARモデ ルを求め
# それぞれの次数に対応した情報量基準を算出する
select.lag <- VARselect(data, type = "const", lag.max = 10)
(nlag <- select.lag$selection[1]) #AIC によるラグ選択

# VAR(y, p = 1, type = c("const", "trend", "both", "none"),
#		season = NULL, exogen = NULL, lag.max = NULL,
#		ic = c("AIC", "HQ", "SC", "FPE"))
# 
# VAR(p)(ベクトル自己回帰モデル)の生成
var.model <- VAR(data, p = nlag) 
eps <- residuals(var.model) # VAR(p=nlag) の残差

# zoo オブジェクトの生成
d.index <- index(data)
eps <- zoo(eps, d.index[(nlag + 1):nobs])
head(eps)

# VARモデルの残差に対して
# Dynamic Conditional Correlation GARCH モデルの当てはめ

# 条件付分散式パラメータ設定
a <- c(0.003, 0.005, 0.001, 0.002, 0.001, 0.003, 0.004, 0.002)
A <- diag(c(0.2, 0.3, 0.15, 0.1, 0.2, 0.2, 0.1, 0.2))
B <- diag(c(0.75, 0.6, 0.8, 0.7, 0.4, 0.6, 0.7, 0.6))
# 条件付相関係数式 パラメータ設定
dcc.para <- c(0.1, 0.8)

# dcc.estimation(inia, iniA, iniB, ini.dcc, dvar, model,
#		method="BFGS", gradient=1, message=1)
#
# DCC GARCH での推定
# 1. 条件付分散式のパラメータについて尤度関数最大化
# 2. 条件付相関関係式のパラメータに関する尤度最大化
#
# 結果:
#    out    : パラメータ推定値とその標準誤差          
#    loglik : パラメータ推定値で評価した対数尤度関数の値 
#    h      : 推定された条件付分散   
#    DCC    : 推定された動的条件付相関係数                        
#    std.resid : 変数(eps)を条件付標準誤差(hの平方根)で除した標準化残差        
#    first  : 第1段階目(条件付分散式に関する対数尤度関数の最適化)の結果
#    second : 第1段階目(条件付相関係数に関する対数尤度関数の最適化)の結果

dcc.results <- dcc.estimation(inia = a, iniA = A, iniB = B, 
				ini.dcc = dcc.para, dvar = eps, model = "diagonal")

dcc.results$out #パラメータ推定値と標準誤差
dcc.results$loglik #パラメータ推定値で評価した対数尤度関数の値 
dcc.results$first #第1段階目(条件付分散式に関する対数尤度関数の最適化)の結果
dcc.results$second #第1段階目(条件付相関係数に関する対数尤度関数の最適化)の結果

dcc.results$h #推定された条件付分散 
dcc.results$DCC #推定された動的条件付相関係数
dcc.results$std.resid #変数(eps)を条件付標準誤差(hの平方根)で除した標準化残差 

# 動的相関係数の出力
# 動的相関係数dcc.results$DCCは、行ごとにt期の条件付相関係数行列を
# ベクトル化した値として保存されているため、
# 全体として(サンプル数)×(次元の2乗)の行列になっている。
head(dcc.results$DCC,3)

# プロットする際には対角成分や重複した列を除外する必要があるため、
# ndim 次 元 の DCC-GARCHモデルにおいては、diam(ndim) (単位行列)を設定し、
# lower.tri関数によって正方行列の対角成分を含まない下三角行列の位置を指定する。
Id <- diag(ndim) #単位行列
head(dcc.results$DCC[, lower.tri(Id)],3) #下三角行列

# zooオブジェクト生成
DCC <- zoo(dcc.results$DCC[, lower.tri(Id)], d.index[(nlag + 1):nobs])
#列名を付ける
cn <- colnames(data)
cnames <- {}
for (i in 1:(length(cn)-1)) {
	for (j in (i+1):length(cn)) {
		cnames <- c(cnames, paste(cn[i], cn[j], sep = "-"))				
	}
}
colnames(DCC) <- cnames

# データ確認 
head(DCC, 3)

# Dynamic Conditional Correlation Estimates表示(一部)
plot(DCC[,1:6], main = "動的条件付相関係数の推定結果")

動的条件付相関係数の推定結果(一部) 表示

f:id:hamadakoichi:20110428034150j:image:w550


参考文献

Rによる時系列分析入門

Rによる時系列分析入門


参考資料

R言語での時系列分析

R言語での Conditional Correlation GARCH Model (ccgarch)

R言語でのVector Autoregressive Model (vars)

トラックバック - http://d.hatena.ne.jp/hamadakoichi/20110427/p1