簡単な回帰モデル plus GARCH(1,1)
options(warn=1,digits=5,messages=1,tidy=TRUE)
suppressPackageStartupMessages(library(Ecdat))
data(package="Ecdat")
# EcdatはR-Projectのサイトから簡単に入手できるので
#このライブラリーに含まれているMoneyUSを使って
#線形回帰式と誤差項がGARCH(1,1)に従うモデルの推定を
#行う。パッケージはrugarchを使う。
#モデルのスペックを指定するときに、外生変数としての
#説明変数をマトリックスの型で与える必要がある。
#ここで推定するモデルは計算例を示すための例示であり
#経済理論的な意味や根拠は全くない。
#推定しようとするモデルは
#cprの差分＝c0＋c1・inflの差分＋誤差項e(t)
#外部の説明変数としてinflの差分が加わり
#誤差項e(t)はGARCH(1,1)に従うと仮定する。
# diff(cpr(t))=c0+c1・diff(infl)+e(t)
#sigma(t)^2=omega+a・e(t-1)^2+b・sigma(t-1)^2
#最尤法でc0,c1,omega,a,bを推定しようとする。
#
dim(MoneyUS)## [1] 164 5# データフレームとしてDDにデータを保存
DD<-as.data.frame(MoneyUS)
head(DD)## m infl cpr y tbr
## 1 1.95138 6.3546 7.4151 2.0367 1.08367
## 2 0.90592 6.3537 7.4136 1.6333 0.81433
## 3 0.35966 6.3626 7.4252 1.3633 0.86967
## 4 3.85845 6.3638 7.4374 1.3100 1.03633
## 5 3.09707 6.3676 7.4631 1.6133 1.25633
## 6 4.69576 6.3615 7.4723 1.9667 1.51433#infl,cprの対前期との差分データを作成
dinf<-diff(DD$infl)
dcpr<-diff(DD$cpr)
length(dcpr)## [1] 163length(dinf)## [1] 163length(DD$infl)## [1] 164#外生変数はマトリックス型で指示する。それ以外の型では
#うまく機能しない。
extdt <- cbind(dinf)
class(extdt)## [1] "matrix"head(extdt)## dinf
## [1,] -0.00093549
## [2,] 0.00896044
## [3,] 0.00118070
## [4,] 0.00381317
## [5,] -0.00610604
## [6,] -0.00493606suppressPackageStartupMessages(library(rugarch))
regspec.garch <- ugarchspec(
variance.model = list(garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0,0),
include.mean = TRUE,
external.regressors = extdt)
)
regspec.garch##
## *---------------------------------*
## * GARCH Model Spec *
## *---------------------------------*
##
## Conditional Variance Dynamics
## ------------------------------------
## GARCH Model : sGARCH(1,1)
## Variance Targeting : FALSE
##
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model : ARFIMA(0,0,0)
## Include Mean : TRUE
## GARCH-in-Mean : FALSE
## Exogenous Regressor Dimension: 1
##
## Conditional Distribution
## ------------------------------------
## Distribution : norm
## Includes Skew : FALSE
## Includes Shape : FALSE
## Includes Lambda : FALSEregspec.garch.fit <- ugarchfit(data = dcpr, spec =
regspec.garch,
out.sample=5)
regspec.garch.fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.006648 0.000760 8.7455 0.000000
## mxreg1 0.173835 0.062814 2.7675 0.005649
## omega 0.000008 0.000000 45.1602 0.000000
## alpha1 0.046390 0.010238 4.5310 0.000006
## beta1 0.855712 0.027047 31.6381 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.006648 0.001051 6.3227 0.000000
## mxreg1 0.173835 0.094377 1.8419 0.065487
## omega 0.000008 0.000000 36.8314 0.000000
## alpha1 0.046390 0.009458 4.9046 0.000001
## beta1 0.855712 0.030649 27.9201 0.000000
##
## LogLikelihood : 519.82
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.5167
## Bayes -6.4197
## Shibata -6.5186
## Hannan-Quinn -6.4773
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 11.90 0.0005609
## Lag[2*(p+q)+(p+q)-1][2] 13.56 0.0002197
## Lag[4*(p+q)+(p+q)-1][5] 14.73 0.0005267
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.0002635 0.9870
## Lag[2*(p+q)+(p+q)-1][5] 1.0096949 0.8575
## Lag[4*(p+q)+(p+q)-1][9] 3.7778561 0.6267
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.2019 0.500 2.000 0.6532
## ARCH Lag[5] 2.2281 1.440 1.667 0.4230
## ARCH Lag[7] 2.8867 2.315 1.543 0.5352
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 18.93
## Individual Statistics:
## mu 0.88570
## mxreg1 0.81192
## omega 2.74758
## alpha1 0.08777
## beta1 0.18683
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.5858 0.5589
## Negative Sign Bias 0.5664 0.5720
## Positive Sign Bias 1.0471 0.2967
## Joint Effect 1.6931 0.6385
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 5.797 0.9984
## 2 30 12.127 0.9975
## 3 40 32.886 0.7439
## 4 50 34.405 0.9434
##
##
## Elapsed time : 0.14981#diff(cpr(t))=0.0066+0.1738・diff(infl)
#sigma(t)^2=0.000+0.046・e(t-1)^2+0.855・sigma(t-1)^2
#とc0,c1,omega,a,bがそれぞれ推定されている。
#回帰モデルの外生変数の係数はmxreg1で示されている。
fcpr <- regspec.garch.fit@fit$fitted.values
plot.ts(fcpr)
sigma2cpr <- regspec.garch.fit@fit$sigma^2
plot.ts(sigma2cpr)
hinf<-diff(DD$infl)
hcpr<-diff(DD$cpr)
length(hinf)## [1] 163length(hcpr)## [1] 163length(fcpr)## [1] 158matplot(cbind(hcpr[1:158],fcpr,hinf[1:158]),type="l",col=c(2:4))
#予測モデル
fmodel <- ugarchforecast(regspec.garch.fit,data=dcpr,n.roll=5,
external.forecasts=list(mregfor=extdt),n.ahead=5)
plot(fmodel, which =c(1))
plot(fmodel, which =c(2))
plot(fmodel, which =c(3))
plot(fmodel, which =c(4))
