chap4_2
qaj
2023-12-09
例子：美元对欧元汇率日对数收益率
- 画出对数收益率平方序列的PACF:
library(readr)
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(xts)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniond.useu <- read_table(
"D:/齐安静 教学/时间序列分析/北大/ftsdata/d-useu9910.txt",
col_types=cols(.default=col_double())
)
xts.useu <- with(d.useu, xts(rate, make_date(year, mon, day)))
xts.useu.lnrtn <- diff(log(xts.useu))[-1,]
eu <- c(coredata(xts.useu.lnrtn))pacf(eu^2, main="")
做图演示
- 对标准化t分布，取\(v=5,10,30\)，作密度图：
dtstd <- function(x, df)
sqrt(df/(df-2)) * dt(x*sqrt(df/(df-2)), df)curve(dtstd(x, 5), -3, 3, lwd=1, lty=1, col="black",
xlab="x", ylab="密度")
curve(dtstd(x, 10), -3, 3, lwd=1, lty=2, col="red", add=TRUE)
curve(dtstd(x, 30), -3, 3, lwd=1, lty=3, col="blue", add=TRUE)
legend("topleft", lty=c(1,2,3), col=c("black", "red", "blue"),
legend=c(expression(v==5), expression(v==10), expression(v==30)))
- 有偏t分布（\(v=5, \xi=0.75, 1,
1.5\)）
dtskew <- function(x, df, xi){
y <- numeric(length(x))
par1 <- gamma((df-1)/2)*sqrt(df-2)/sqrt(pi)/gamma(df/2)*(xi - 1/xi)
par2 <- sqrt(xi^2 + 1/xi^2 - 1 - par1^2)
c1 <- 2*par2/(xi + 1/xi)
sele <- x < -par1/par2
y[sele] <- c1 * dtstd(xi*(par2*x[sele] + par1), df)
y[!sele] <- c1 * dtstd((par2*x[!sele] + par1)/xi, df)
y
}curve(dtskew(x, 5, 1.0), -3, 3, ylim=c(0, 0.6),
lwd=1, lty=1, col="black",
xlab="x", ylab="密度")
curve(dtskew(x, 5, 0.75), -3, 3, lwd=1, lty=2, col="red", add=TRUE)
curve(dtskew(x, 5, 1.5), -3, 3, lwd=1, lty=3, col="blue", add=TRUE)
legend("topleft", lty=c(1,2,3), col=c("black", "red", "blue"),
legend=c(expression(xi==1), expression(xi==0.75), expression(xi==1.5)))
dtskew <- function(x, df, xi){
y <- numeric(length(x))
par1 <- gamma((df-1)/2)*sqrt(df-2)/sqrt(pi)/gamma(df/2)*(xi - 1/xi)
par2 <- sqrt(xi^2 + 1/xi^2 - 1 - par1^2)
c1 <- 2*par2/(xi + 1/xi)
sele <- x < -par1/par2
y[sele] <- c1 * dtstd(xi*(par2*x[sele] + par1), df)
y[!sele] <- c1 * dtstd((par2*x[!sele] + par1)/xi, df)
y
}curve(dtskew(x, 5, 1.0), -3, 3, ylim=c(0, 0.6),
lwd=1, lty=1, col="black",
xlab="x", ylab="密度")
curve(dtskew(x, 5, 0.75), -3, 3, lwd=1, lty=2, col="red", add=TRUE)
curve(dtskew(x, 5, 1.5), -3, 3, lwd=1, lty=3, col="blue", add=TRUE)
legend("topleft", lty=c(1,2,3), col=c("black", "red", "blue"),
legend=c(expression(xi==1), expression(xi==0.75), expression(xi==1.5)))
Intel公司股票ARCH建模实例
- 继续使用1973年到2009年Intel公司股票的月对数收益率。 ARCH效应检验一节进行了ARCH效应检验， 证明有ARCH效应。
- 数据读入
d.intel <- read_table("D:/齐安静 教学/时间序列分析/北大/ftsdata/m-intcsp7309.txt", col_types=cols(.default=col_double(),
date=col_date("%Y%m%d")))
xts.intel <- xts(log(1 + d.intel[["intc"]]), d.intel[["date"]])
tclass(xts.intel) <- "yearmon"
ts.intel <- ts(c(coredata(xts.intel)), start=c(1973,1), frequency=12)
at <- ts.intel - mean(ts.intel)forecast::Acf(at^2, lag.max=36, main="")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
- 残差的平方的PACF:
pacf(at^2, lag.max=36, main="")
- 使用fGarch包的
garchFit()函数建立ARCH模型。
library(fGarch)## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")mod1 <- garchFit(~ 1 + garch(3,0), data=c(ts.intel), trace=FALSE)summary(mod1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(3, 0), data = c(ts.intel), trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(3, 0)
## <environment: 0x00000222ac8e0670>
## [data = c(ts.intel)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3
## 0.012567 0.010421 0.232889 0.075069 0.051994
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.012567 0.005515 2.279 0.0227 *
## omega 0.010421 0.001238 8.418 <2e-16 ***
## alpha1 0.232889 0.111541 2.088 0.0368 *
## alpha2 0.075069 0.047305 1.587 0.1125
## alpha3 0.051994 0.045139 1.152 0.2494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 303.9607 normalized: 0.6845963
##
## Description:
## Thu Dec 14 15:26:30 2023 by user: qaj
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 203.3620090 0.000000e+00
## Shapiro-Wilk Test R W 0.9635971 4.898648e-09
## Ljung-Box Test R Q(10) 9.2607810 5.075464e-01
## Ljung-Box Test R Q(15) 19.3674758 1.975620e-01
## Ljung-Box Test R Q(20) 20.4698241 4.289060e-01
## Ljung-Box Test R^2 Q(10) 7.3221345 6.947236e-01
## Ljung-Box Test R^2 Q(15) 27.4153158 2.552912e-02
## Ljung-Box Test R^2 Q(20) 28.1511237 1.058699e-01
## LM Arch Test R TR^2 25.2334658 1.375448e-02
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.346670 -1.300546 -1.346920 -1.328481- 因为结果中\(\alpha_2\)和\(\alpha_3\)的估计值是不显著的， 可拟合ARCH(1)模型为波动率方程：
mod2 <- garchFit( ~ 1 + garch(1,0), data=c(ts.intel), trace=FALSE)summary(mod2)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(1, 0), data = c(ts.intel), trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(1, 0)
## <environment: 0x00000222aebe5f90>
## [data = c(ts.intel)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1
## 0.013130 0.011046 0.374976
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.013130 0.005318 2.469 0.01355 *
## omega 0.011046 0.001196 9.238 < 2e-16 ***
## alpha1 0.374976 0.112620 3.330 0.00087 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 299.9247 normalized: 0.675506
##
## Description:
## Thu Dec 14 15:26:30 2023 by user: qaj
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 144.3781954 0.000000e+00
## Shapiro-Wilk Test R W 0.9678175 2.670339e-08
## Ljung-Box Test R Q(10) 12.1224742 2.769434e-01
## Ljung-Box Test R Q(15) 22.3070432 1.000021e-01
## Ljung-Box Test R Q(20) 24.3341140 2.281020e-01
## Ljung-Box Test R^2 Q(10) 16.5780431 8.423799e-02
## Ljung-Box Test R^2 Q(15) 37.4434390 1.089753e-03
## Ljung-Box Test R^2 Q(20) 38.8138993 7.031666e-03
## LM Arch Test R TR^2 27.3289440 6.926884e-03
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.337499 -1.309824 -1.337589 -1.326585- 计算标准化残差
resi <- residuals(mod2, standardize=TRUE)- 标准化残差的时序图：
plot(mod2, which=9)
- 标准化残差\(\{ \tilde a_t \}\)的ACF:
plot(mod2, which=10)
- \(\{\tilde a_t^2\}\)的ACF:
plot(mod2, which=11)
- \(\{\tilde a_t^2\}\)的PACF:
pacf(resi^2, lag.max=36, main="")
- 标准化残差相对于标准正态分布的QQ图
plot(mod2, which=13)
- 显示拟合的条件标准差序列（波动率序列）：
plot(mod2, which=2)
- 获得拟合的波动率用volatility()函数：
mod2v <- volatility(mod2)
head(mod2v)## [1] 0.1306624 0.1051189 0.1450052 0.1101684 0.1134644 0.1294741- 获得拟合的均值用
fitted()函数。 按理说这个模型的均值方程应该是常数0.0131
mod2f <- fitted(mod2)
head(mod2f)## 1 2 3 4 5 6
## 0.01313002 0.01313002 0.01313002 0.01313002 0.01313002 0.01313002- 这里返回的是原始\(r_t\)的序列。
head(ts.intel)## Jan Feb Mar Apr May
## 1973 0.009999835 -0.150012753 0.067064079 0.082948635 -0.110348491
## Jun
## 1973 0.125162849- 用
predict()函数作超前若干步的预测，如：
mod2p <- predict(mod2, n.ahead=12)
mod2p## meanForecast meanError standardDeviation
## 1 0.01313002 0.1090512 0.1090512
## 2 0.01313002 0.1245214 0.1245214
## 3 0.01313002 0.1298481 0.1298481
## 4 0.01313002 0.1317900 0.1317900
## 5 0.01313002 0.1325108 0.1325108
## 6 0.01313002 0.1327801 0.1327801
## 7 0.01313002 0.1328809 0.1328809
## 8 0.01313002 0.1329187 0.1329187
## 9 0.01313002 0.1329329 0.1329329
## 10 0.01313002 0.1329382 0.1329382
## 11 0.01313002 0.1329402 0.1329402
## 12 0.01313002 0.1329410 0.1329410- 原始序列的均值
mean(ts.intel)## [1] 0.0143273- 原始序列的标准差
sd(ts.intel)## [1] 0.1269101Intel股票问题改用t分布
- 在fGarch的
garchFit()函数中加选项cond.dist="std"。
mod3 <- garchFit(~ 1 + garch(1,0), data=ts.intel,
cond.dist="std", trace=FALSE)- 结果
summary(mod3)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(1, 0), data = ts.intel, cond.dist = "std",
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(1, 0)
## <environment: 0x00000222a8b6a308>
## [data = ts.intel]
##
## Conditional Distribution:
## std
##
## Coefficient(s):
## mu omega alpha1 shape
## 0.017202 0.011816 0.277476 5.970256
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.017202 0.005195 3.311 0.000929 ***
## omega 0.011816 0.001560 7.574 3.62e-14 ***
## alpha1 0.277476 0.107183 2.589 0.009631 **
## shape 5.970256 1.529522 3.903 9.49e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 315.0899 normalized: 0.709662
##
## Description:
## Thu Dec 14 15:26:31 2023 by user: qaj
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 157.7799623 0.000000e+00
## Shapiro-Wilk Test R W 0.9663975 1.488220e-08
## Ljung-Box Test R Q(10) 12.8594024 2.316394e-01
## Ljung-Box Test R Q(15) 23.4063309 7.588549e-02
## Ljung-Box Test R Q(20) 25.3740038 1.874954e-01
## Ljung-Box Test R^2 Q(10) 19.9609448 2.962426e-02
## Ljung-Box Test R^2 Q(15) 42.5555179 1.845072e-04
## Ljung-Box Test R^2 Q(20) 44.0674224 1.473957e-03
## LM Arch Test R TR^2 29.7607256 3.033495e-03
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.401306 -1.364407 -1.401466 -1.386755- 作标准化残差相对于标准化t分布的QQ图：
plot(mod3, which=13)
欧元汇率ARCH建模实例
- 考虑1999-01-04到2010-08-20的欧元对美元汇率的日对数收益率。显示数据有ARCH效应。 该序列是纯条件异方差模型的典型例子。
d.useu <- read_table(
"D:/齐安静 教学/时间序列分析/北大/ftsdata/d-useu9910.txt",
col_types=cols(.default=col_double())
)
xts.useu <- with(d.useu, xts(rate, make_date(year, mon, day)))
xts.useu.lnrtn <- diff(log(xts.useu))[-1,]
eu <- c(coredata(xts.useu.lnrtn))- \(a_t^2\)的ACF
forecast::Acf(eu^2, lag.max=36, main="")
- \(a_t^2\)的PACF:
pacf(eu^2, lag.max=36, main="")
- 在低阶到滞后11可以， 但是高阶仍有较大的值。 考虑建立ARCH(11)模型。 采用条件高斯似然函数：
mod4 <- garchFit(~ 1 + garch(11,0), data=eu, trace=FALSE)summary(mod4)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(11, 0), data = eu, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(11, 0)
## <environment: 0x00000222ae09e8e0>
## [data = eu]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3 alpha4
## 1.2645e-04 1.8903e-05 1.6607e-02 4.4563e-02 2.7212e-02 8.0370e-02
## alpha5 alpha6 alpha7 alpha8 alpha9 alpha10
## 5.0111e-02 9.2190e-02 7.5282e-02 6.9537e-02 3.3466e-02 2.7823e-02
## alpha11
## 3.8773e-02
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 1.265e-04 1.110e-04 1.140 0.254446
## omega 1.890e-05 1.727e-06 10.945 < 2e-16 ***
## alpha1 1.661e-02 1.575e-02 1.055 0.291577
## alpha2 4.456e-02 2.085e-02 2.137 0.032593 *
## alpha3 2.721e-02 1.700e-02 1.601 0.109353
## alpha4 8.037e-02 2.362e-02 3.402 0.000669 ***
## alpha5 5.011e-02 2.127e-02 2.355 0.018500 *
## alpha6 9.219e-02 2.274e-02 4.053 5.05e-05 ***
## alpha7 7.528e-02 2.406e-02 3.129 0.001755 **
## alpha8 6.954e-02 2.455e-02 2.832 0.004622 **
## alpha9 3.347e-02 2.022e-02 1.655 0.097835 .
## alpha10 2.782e-02 1.820e-02 1.528 0.126412
## alpha11 3.877e-02 1.906e-02 2.035 0.041894 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 10698.47 normalized: 3.652603
##
## Description:
## Thu Dec 14 15:26:45 2023 by user: qaj
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 360.8032273 0.000000e+00
## Shapiro-Wilk Test R W 0.9891735 3.894032e-14
## Ljung-Box Test R Q(10) 15.7762969 1.062177e-01
## Ljung-Box Test R Q(15) 21.5010500 1.215699e-01
## Ljung-Box Test R Q(20) 24.7744465 2.101968e-01
## Ljung-Box Test R^2 Q(10) 4.8009795 9.040700e-01
## Ljung-Box Test R^2 Q(15) 16.7306668 3.352185e-01
## Ljung-Box Test R^2 Q(20) 27.5606043 1.202156e-01
## LM Arch Test R TR^2 11.9678220 4.482674e-01
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -7.296329 -7.269777 -7.296368 -7.286766