Rソースコード

GARCH(1,1)のデータ生成とパラメータ推定

GARCH(1,1)の時系列を人工的に作り図示してみよう. また, この時系列のパラメータ推定を行ってみよう. ここでは R-Forge(http://r-forge.r-project.org/R/?group_id=156)よりGARCHモデル推定ライブラリfGarchをダウンロードして用いた. fGarchにはGARCH時系列生成関数, garchSim(spec, n, n.stat, extended), パラメータ推定関数 garchFit(formula, data)などが含まれている.

GARCH(1,1)のデータ生成とパラメータ推定


library(fGarch)
library(FinTS)
omega = 1e-06
alpha = 0.1
beta = 0.8
spec = garchSpec(model = list(omega, alpha, beta))
x <- garchSim(spec, n = 1000)
ts.plot(x)
res <- garchFit( formula = ~garch(1, 1), data = x )
cat(sprintf("%f %f %f\n",omega,alpha,beta))
cat(sprintf("%s\n",attributes(res)$fit$coef))

インデックスに戻る


ACD(1,1)のデータ生成とパラメータ推定

ACD(1,1)の時系列を人工的に作り図示し, この時系列のパラメータ推定を行ってみよう. ここでは, R-Forge(http://r-forge.r-project.org/R/?group_id=156) よりACD推定ライブラリfACDをダウンロードして用いた. fACDにはACD時系列生成関数, ACD_Simul(nr, Coeff, distrib, typeACD), パラメータ推定関数 ACD_Fit(x, qLag, pLag, distrib, typeACD)などが含まれている (本書執筆現在, fACDはR3.0.0以降のバージョンではサポートされなくなった関数を含むためそれ以前のバージョンでのみで作動することに注意されたい).

EACD(1,1)のデータ生成とパラメータ推定


library(fACD)
library(FinTS)
C <- array(1:1);
C$w = 0.1
C$p = 0.5
C$q = 0.3
x <- ACD_Simul(1000, C, "exp", "ACD")
ts.plot(x)
res <- ACD_Fit(x,1,1,"exp","ACD")
CC <- attributes(res)$Coeff
AIC = -2.0*attributes(res)$LL
AIC = AIC +2.0*attributes(res)$sizeModel$nParameters
cat(sprintf("%f %f %f %f\n",CC$w,CC$p,CC$q,AIC))

インデックスに戻る


yuimaパッケージ

ここでは, 三菱UFJフィナンシャルグループ(8306) とみずほフィナンシャルグループ(8411) の2つの株式の日次相関係数(T=1)を, 日経NEEDSデータを使って推定してみよう. 分析データは約定価格系列, データ期間は, 2010年1年間, 時間解像度は秒である. 分析の前に, 同一のタイムスタンプを持つレコードは, 最初のもののみを残し, 2番目以降は削除した. 推定法としては, HY推定量をマイクロストラクチャ・ノイズのある価格系列に対応できるように改良したPre-averaged Hayashi-Yoshida推定量("PHY", (Christensen et al., 2010)), 上記 (Kunitomo and Sato, 2013)によるSIML法("SIML"), 同じく実現カーネル法("RK", (Barndorff-Nielsen et al., 2011)), さらに, 上記 (Zhang et al., 2005) による実現ボラティリティ推定法を多変量に拡張した2スケール共分散推定法 (Zhang, 2011)の4つを同時に使用し, その大きさを比較する. Ryuimaパッケージに含まれるcce関数を用いれば, これらは簡単に計算することができる. (yuima開発プロジェクトは, 東京大学大学院数理科学研究科の吉田朋広教授のグループによって行われている.)

yuimaパッケージによる銘柄間の相関係数の推定


library(yuima)
# define yuima calss data
# x.logprcは対数価格ベクトル, x.hmsはタイムスタンプ
xprc <- zoo(x.logprc,order.by=x.hms)	
# zoo()によって, 非等間隔時系列データを作成
yprc <- zoo(y.logprc,order.by=y.hms) 	
# 2つの非等間隔時系列データを結合(マージ) 
xyprc<- merge.zoo(xprc,yprc)		
# setData()によって, yuimaクラスとして登録
xyprc.yuima <- setData(xyprc) 		
# calculate correlation
# PHY法,  タイムスタンプは1秒単位(utime=1)
cce(xyprc.yuima,method="PHY",utime=1)	
cce(xyprc.yuima,method="SIML",utime=1) 	# SIML法
cce(xyprc.yuima,method="RK",c.RK=1, utime=1)	# RK法
cce(xyprc.yuima,method="TSCV",c.two=1, utime=1)	# TSCV法

 [文献]
Christensen, K., Kinnebrock, S. and Podolskij, M. (2010) "Pre-averaging estimators of the ex-post covariance matrix in noisy diffusion models with non-synchronous data", J. Econometrics, Vol. 159, pp. 116-133.
Kunitomo, N. and Sato, S. (2013) "Separating Information Maximum Likelihood estimation of the integrated volatility and covariance with micro-market noise", North Amer. J. Econ. Fin..
Barndorff-Nielsen, O. E., Hansen, P. R., Lunde, A. and Shephard, N. (2011) "Multivariate realised kernels: consistent positive semi-definite estimators of the covariation of equity prices with noise and non-synchronous trading", J. Econometrics, Vol. 162(2), pp. 149-169.
Zhang, L., Mykland, P. A. and Ait-Sahalia, Y. (2005) "A Tale of Two Time Scales: Determining Integrated Volatility with Noisy High-frequency Data", J. Amer. Statist. Assoc., Vol. 100, pp. 1394-1411.
Zhang, L. (2011) "Estimating covariation: Epps effect, microstructure noise", J. Econometrics, Vol. 160, pp. 33-47.

インデックスに戻る


正規性の検定

統計言語Rにおいて人工的に作成した300ポイントの正規乱数からKS検定を行うためのサンプルソースコードを以下に示す.

Kolmogorov-Smirnov検定による正規性の確認


pks <- function(x){
  i<-1
  cdf <- 0.0
  while (1) {
    d <- 2.0*exp(-2.0*(2*i-1)*(2*i-1)*x*x)
    d <- d - 2.0*exp(-2.0*(2*i)*(2*i)*x*x)
    if(d < 1e-16) break
    cdf <- cdf + d
    i <- i+1
  }
  cdf
}
KStest <- function(x){
  x <- x[!is.na(x)]
  n <- length(x)
  PVAL <- NULL
  sx <- sort(x)
  xx <- c(1:n)
  for(i in 1:n){
    xx[i] <- length(sx[sx<sx[i]])/length(x)
  }
  y <- pnorm(sx,mean(sx),sd(sx))
  D <- max(abs(y-xx)*sqrt(n))
  PVAL <- pks(D) 
  cat(sprintf("D=%f, p=%f\n",D,PVAL))
}
x<-rnorm(300)
KStest(x)

α=0.05 (5%有意水準)とすると, 得られるp値が0.05以下である場合に帰無仮説H0 (2つの分布は一致する) は棄却され, 標本データが従う分布は正規分布ではないと結論付けることができる (コードを実行して読者自身で値を確かめてほしい).

インデックスに戻る


(C)2016 Takaki Hayashi and Aki-Hiro Sato. All rights reserved. 無断転用を禁ず