chap2_5

案例：3M公司股票月对数收益率（1946-2008）

library(readr)
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(xts)
d <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/m-3m4608.txt",
  col_types=cols(.default=col_double(),
                 date=col_date(format="%Y%m%d")))
mmm <- xts(log(1 + d[["rtn"]]), d$date)
rm(d)
tclass(mmm) <- "yearmon"
ts.3m <- ts(coredata(mmm), start=c(1946,2), frequency=12)
ts.3mlogp <- ts(cumsum(c(c(ts.3m))), start=c(1946,2), frequency=12)
c(mean=mean(ts.3m), sd=sd(ts.3m))

##       mean         sd 
## 0.01029941 0.06371910

对数收益率的白噪声检验:

Box.test(ts.3m, type="Ljung", lag=12)

## 
##  Box-Ljung test
## 
## data:  ts.3m
## X-squared = 27.688, df = 12, p-value = 0.006143

\(\{p_t\}\)的ACF图形如下：

forecast::Acf(ts.3mlogp, main="", lag=36)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

图形分解

plot(c(time(ts.3mlogp)), c(ts.3mlogp), ylim=c(-1, 8.2), type="l", 
     xlab="Year", ylab="ln(Price)")
tmp.x <- seq(length(ts.3mlogp))
tmp.y <- 0.01030*tmp.x
tmp.y2 <- c(ts.3mlogp) - tmp.y
tmp.y <- ts(tmp.y, start=c(1946,2), frequency=12)
lines(c(time(tmp.y)), c(tmp.y), col="red", lwd=2, lty=3)
lines(c(time(tmp.y)), c(tmp.y2), col="green")
legend("topleft", lty=c(1,3,1), col=c("black", "red", "green"),
       legend=c("log(P)", "Linear Trend", "Random Walk"))
abline(h=0, col="gray", lty=3)

删除临时变量

rm(tmp.x, tmp.y, tmp.y2)

案例2: 美国的国民生产总值(GNP)

读入数据，对数化处理，并做一阶差分

da <- read_table("D:/齐安静 教学/时间序列分析/北大/ftsdata/q-gnp4710.txt", col_types=cols(
  .default = col_double()))
gnp <- ts(log(da[["VALUE"]]), 
          start=c(1947, 1), frequency=4)
dgnp <- diff(gnp)
rm(da)

GNP的对数值的图形：

plot(gnp, xlab="year", ylab="log(GNP)")

对数GNP明显地不平稳。其ACF图形如下：

forecast::Acf(gnp, main="")

对数GNP差分的时序图

plot(dgnp, xlab="year", ylab="log(GNP)")
abline(h=0, col="gray", lty=3)

对数GNP差分的ACF

forecast::Acf(dgnp, main="")

对数GNP差分的PACF

forecast::Pacf(dgnp, main="")

如果用AR(\(p\))为对数GNP差分序列建模， PACF图提示用\(p=16\)。用AIC定阶：

ar(dgnp, method="mle")

## 
## Call:
## ar(x = dgnp, method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.4318   0.1985  -0.1180   0.0189  -0.1607   0.0900   0.0615  -0.0814  
##       9  
##  0.1940  
## 
## Order selected 9  sigma^2 estimated as  8.918e-05

ADF的基础模型需要一个AR阶数，取\(p=9\)。用fUnitRoots::adfTest()对GNP的对数值进行ADF单位根检验：

fUnitRoots::adfTest(gnp, lags=9, type="c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 9
##   STATISTIC:
##     Dickey-Fuller: -1.8467
##   P VALUE:
##     0.3691 
## 
## Description:
##  Tue Nov  7 13:17:56 2023 by user: qaj

结果p值较大，说明不能拒绝零假设，即对数GNP序列有单位根。
在fUnitRoots::adfTest()中，用lag=指定检验所用AR的阶数，用type="c"指定基础模型允许有一个非零常数项，用type="nc"指定基础模型不允许有任何的常数项和线性项，用type="ct"指定基础模型允许有常数项和线性项。
GNP对数序列的图形也像是有非随机线性增长趋势的情况。为此，仍使用ADF检验，但是允许有非随机常数项和线性项：

fUnitRoots::adfTest(gnp, lags=9, type="ct")

## Warning in fUnitRoots::adfTest(gnp, lags = 9, type = "ct"): p-value greater
## than printed p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 9
##   STATISTIC:
##     Dickey-Fuller: -0.0094
##   P VALUE:
##     0.99 
## 
## Description:
##  Tue Nov  7 13:17:56 2023 by user: qaj

尝试人为地拟合非随机线性增长趋势，检验残差是否有单位根：

tmp.t <- c(time(gnp))
tmp.y <- residuals( lm(c(gnp) ~ tmp.t) )
fUnitRoots::adfTest(tmp.y, type="nc")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.0763
##   P VALUE:
##     0.592 
## 
## Description:
##  Tue Nov  7 13:17:56 2023 by user: qaj

使用urca包的ur.df()执行单位根ADF检验，基础模型带有常数漂移项

library(urca)

## Warning: package 'urca' was built under R version 4.3.1

summary(ur.df(gnp, type="drift", lags=9))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.033589 -0.004595  0.000503  0.003911  0.035557 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0135157  0.0045219   2.989  0.00310 ** 
## z.lag.1     -0.0009040  0.0004895  -1.847  0.06606 .  
## z.diff.lag1  0.3903745  0.0641690   6.084 4.81e-09 ***
## z.diff.lag2  0.2057039  0.0681595   3.018  0.00283 ** 
## z.diff.lag3 -0.1111342  0.0693754  -1.602  0.11053    
## z.diff.lag4  0.0263875  0.0692327   0.381  0.70345    
## z.diff.lag5 -0.1399063  0.0684739  -2.043  0.04216 *  
## z.diff.lag6  0.0857675  0.0704276   1.218  0.22453    
## z.diff.lag7  0.0520310  0.0696968   0.747  0.45610    
## z.diff.lag8 -0.0879122  0.0685969  -1.282  0.20127    
## z.diff.lag9  0.1816295  0.0639332   2.841  0.00490 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.009396 on 232 degrees of freedom
## Multiple R-squared:  0.3142, Adjusted R-squared:  0.2847 
## F-statistic: 10.63 on 10 and 232 DF,  p-value: 8.56e-15
## 
## 
## Value of test-statistic is: -1.8467 6.7291 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

可以将统计量的第一个(-1.8467)与tau2统计量的5pct临界值(-2.87)比较，没有小于临界值，所以可以接受零假设，认为有单位根。
使用tseries包的adf.test()执行单位根ADF检验：

tseries::adf.test(gnp, k=9)

## Warning in tseries::adf.test(gnp, k = 9): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  gnp
## Dickey-Fuller = -0.0093764, Lag order = 9, p-value = 0.99
## alternative hypothesis: stationary

案例：标普500（1950-2008）

读入数据

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

da <- read_table("D:/齐安静 教学/时间序列分析/北大/ftsdata/d-sp55008.txt", col_types=cols(.default = col_double()))
xts.sp5d <- xts(da[,-(1:3)], make_date(da$year, da$mon, da$day))
str(xts.sp5d)

## An xts object on 1950-01-03 / 2008-04-11 containing: 
##   Data:    double [14662, 6]
##   Columns: open, high, low, close, volume ... with 1 more column
##   Index:   Date [14662] (TZ: "UTC")

对数处理，一阶差分

sp5d <- log(xts.sp5d[,"close"])
delta.sp5d <- diff(coredata(sp5d)[,1])

日对数收盘价的时间序列图：

plot(sp5d, type="l", main="S&P 500 Daily Log Close",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto")

一阶差分的时间序列图

plot(diff(sp5d), type="l", main="S&P 500 Daily Log Return",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto")

一阶差分的PACF

pacf(delta.sp5d, main="")

用AIC对日对数收益率定阶

ar(delta.sp5d, method="mle")

## 
## Call:
## ar(x = delta.sp5d, method = "mle")
## 
## Coefficients:
##       1        2  
##  0.0721  -0.0387  
## 
## Order selected 2  sigma^2 estimated as  8.068e-05

AIC定阶\(p=2\)。按\(p=2\)对标普500日对数收盘价作ADF白噪声检验：

fUnitRoots::adfTest(c(coredata(sp5d)), lags=2, type="ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.0179
##   P VALUE:
##     0.5708 
## 
## Description:
##  Tue Nov  7 13:17:59 2023 by user: qaj

对日对数收益率（即对数收盘价序列的差分）进行ADF检验

fUnitRoots::adfTest(delta.sp5d, lags=2, type="ct")

## Warning in fUnitRoots::adfTest(delta.sp5d, lags = 2, type = "ct"): p-value
## smaller than printed p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -70.5501
##   P VALUE:
##     0.01 
## 
## Description:
##  Tue Nov  7 13:17:59 2023 by user: qaj

chap2_5

qaj

2023-11-07

案例2: 美国的国民生产总值(GNP)

案例：标普500（1950-2008）