北大金融时间序列分析讲义第3章： 线性时间序列模型

例子：苹果公司2007年到2017年股票日收盘价

chartSeries(
  AAPL, type="line", TA=NULL,
  subset="2003/2017",
  major.ticks="years", minor.ticks=FALSE,
  theme="white", name="Apple"
  )

股价序列呈现缓慢的、非单调的上升趋势，
局部又有短暂的波动。

例子：可口可乐公司盈利季度数据

可口可乐公司每季度发布的每股盈利数据。
读入：

da <- read_table(
  "q-ko-earns8309.txt",
  col_types=cols(
    .default = col_double(),
    pends = col_date("%Y%m%d"),
    anntime = col_date("%Y%m%d")
    ) )
ko.Rqtr <- xts(da[["value"]], da[["pends"]])

时间序列图：

chartSeries(
  ko.Rqtr, type="line", TA=NULL,
  major.ticks="years", minor.ticks=FALSE,
  theme="white", name="Coca Kola Quarterly Return"
  )

序列仍体现出缓慢的、非单调的上升趋势，又有明显的每年的周期变化（称为季节性），
还有短期的波动。

下面用基本R的plot()作图并用不同颜色标出不同季节。

tmp.x <- year(index(ko.Rqtr)) + (quarter(index(ko.Rqtr))-1)/4
tmp.y <- c(coredata(ko.Rqtr))
plot(tmp.x, tmp.y, type="l", col="gray",
     xlab="year", ylab="Return")
cpal <- c("green", "red", "yellow", "black")
points(tmp.x, tmp.y, pch=16, 
       col=cpal[quarter(index(ko.Rqtr))])
legend("topleft", pch=16, col=cpal,
       legend=c("Spring", "Summer", "Autumn", "Winter"))

用ggplot2包绘图：

da_b <- da %>%
  mutate(time = as.yearqtr(pends),
    qtr = factor(quarter(pends),
       levels = 1:4,
       labels=c("Spring", "Summer", "Autumn", "Winter")))
ggplot(data = da_b, mapping = aes(
  x = time, y = value)) +
  geom_point(mapping=aes(color=qtr)) +
  geom_line()

现在可以看出，每年一般冬季和春季最低，
夏季最高，秋季介于夏季和冬季之间。

例子：标普500指数月对数收益率

d <- read_table(
  "m-ibmsp-2611.txt",
  col_types=cols(.default=col_double(),
                 date=col_date(format="%Y%m%d")))
sp.rmon <- xts(log(1 + d[["sp"]]), d[["date"]])

chartSeries(
  sp.rmon, type="line", TA=NULL,
  major.ticks="auto", minor.ticks=FALSE,
  theme="white", name="S&P 500 Monthly Returns"
  )

收益率在0上下波动，除了个别时候基本在某个波动范围之内。

例子：美国国债3月期和6月期周利率

d <- read_table2(
  "w-tb3ms.txt",
  col_types=cols(.default=col_double()))

## Warning: `read_table2()` was deprecated in readr 2.0.0.
## Please use `read_table()` instead.

x1 <- xts(d[["rate"]], make_date(d[["year"]], d[["mon"]], d[["day"]]))
d <- read_table2(
  "w-tb6ms.txt",
  col_types=cols(.default=col_double()))

## Warning: `read_table2()` was deprecated in readr 2.0.0.
## Please use `read_table()` instead.

x2 <- xts(d[["rate"]], make_date(d[["year"]], d[["mon"]], d[["day"]]))
tb36ms <- merge(x1, x2)
names(tb36ms) <- c("tb3ms", "tb6ms")
rm(d, x1, x2)

用xts包的plot()函数作图：

plot(tb36ms, type="l", grid.ticks.on="years")

聚焦到2004年的数据：

plot(tb36ms, type="l", subset="2004")

红色是6月期国债利率，
黑色是3月期国债。
一般6月期高，
但是有些时期3月期超过了6月期，如1980年：

plot(tb36ms, type="l", subset="1980")

相关系数

图3.5是IBM股票月度简单收益率对标普500收益率的散点图。
从图中看出，
两者有明显的正向相关关系。

两个随机变量X和Y的相关系数定义为

ρ(X,Y)=ρxy=Cov(X,Y)Var(X)Var(Y)‾‾‾‾‾‾‾‾‾‾‾‾‾√=E[(X−μx)(Y−μy)]E(X−μx)2E(Y−μy)2‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√

如果有(X,Y)的独立同分布样本(xt,yt)，
t=1,2,…,T，
可估计相关系数为

ρ̂ xy=∑Tt=1(xt−x¯)(yt−y¯)∑Tt=1(xt−x¯)2∑Tt=1(yt−y¯)2‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√

对于不独立的样本，
比如时间序列样本，
也可以计算相关系数，
其估计合理性需要一些模型假设。

对于联合分布非正态的情况，
有时相关系数不能很好地反映X和Y的正向或者负向的相关。
斯皮尔曼（Spearman）相关系数是计算X的样本的秩（名次）与Y的样本的秩之间的相关系数，
也称为Spearman rank correlation。

另一种常用的非参数相关系数是肯德尔tau(Kendall’s τ)系数，
反映了一致数对和非一致数对之间的差别。
对随机向量(X,Y)，
设(X1,Y1), (X2,Y2)相互独立且联合分布与(X,Y)联合分布相同，
定义X和Y的肯德尔tau系数为

τ=P[(X1−X2)(Y1−Y2)>0]−P[(X1−X2)(Y1−Y2)<0]

即两个观测的分量次序一致的概率减去分量次序相反的概率。
一致的概率越大，说明两个的正向相关性越强。

对IBM收益率与标普收益率数据计算这三种相关系数：

## [1] 0.6395979

cor(d[,"sp"], d[,"ibm"], method="spearman")

## [1] 0.6065789

cor(d[,"sp"], d[,"ibm"], method="kendall")

## [1] 0.4328066

自相关函数与白噪声

设{Xt}为弱平稳序列，
{γk}为自协方差函数。
则

ρ(Xt−k,Xt)=Cov(Xt−k,Xt)Var(Xt−k)Var(Xt)‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√=γkγ0γ0‾‾‾‾√=γkγ0, k=0,1,…, ∀t

记ρk=γk/γ0，
这是Xt−k与Xt的相关系数且与t无关，
称{ρk,k=0,1,…}为时间序列{Xt}的自相关函数
（Autocorrelation function, ACF）。
ρ0=1。

如果弱平稳序列{Xt}满足ρk=0, k=1,2,…，
称{Xt}为白噪声序列。
如果随机变量序列{Xt}独立且期望和方差不随时间而变，
则{Xt}是白噪声序列，
称为独立白噪声。
如果独立白噪声还是同分布的，
称为独立同分布白噪声。

适当条件下ρk可以从时间序列样本估计为

ρ̂ k=γ̂ kγ̂ 0, k=0,1,…

ρ̂ 0=1。
称ρ̂ k,k=1,2,…为样本自相关函数。

如果时间序列严平稳遍历，
则ρ̂ k是ρk的强相合估计。

若{Xt}为有二阶矩的独立同分布随机变量列，
则
ρ̂ k(k>0)渐近服从N(0,1T)。

如果{εt}是零均值独立同分布白噪声，
q为非负整数，
{ψj,j=0,1,…,q}是数列，ψ0=1，

Xt=μ+∑j=0qψjεt−j, t∈ℤ,

则从{Xt,t=1,…,T}计算的ACF满足:
当k>q时，
T‾‾√ρ̂ k渐近服从N(0,1+2∑qj=0ρ2j),
这称为Bartlett公式。
参见 (何书元 2003) P.131 §4.2的例2.1。
原始文献：
MAURICE STEVENSON BARTLETT,
On the Theoretical Specification and Sampling Properties of Auto-Correlated Time Series,
Journal of the Royal Statistical Society (Supplement) 8 (1946), pp. 24-41.

在基本R软件中，
acf(x)可以估计时间序列x的自相关函数并对其前面若干项画图。

例3.1 例：CRSP的第10分位组合的月对数收益率，
1967-1到2009-12。
第10分位组合是NYSE、AMEX、NASDAQ市值最小的10%股票组成的投资组合，
每年都重新调整。

• CRSP是Center for Research in Security Prices,
  位于Chicago Booth。
• NYSE(The New York Stock Exchange, 纽约证券交易所),
• AMEX(American Stock Exchange, 美国证券交易所，在纽约华尔街附近)，
• NASDAQ(National Association of Securities Dealers Automated Quotations，纳斯达克，位于纽约)。

d <- read_table2(
  "m-dec12910.txt", col_types=cols(
    .default=col_double(),
    date=col_date(format="%Y%m%d")))

## Warning: `read_table2()` was deprecated in readr 2.0.0.
## Please use `read_table()` instead.

dec <- xts(as.matrix(d[,-1]), d$date)
tclass(dec) <- "yearmon"
d10 <- ts(coredata(dec)[,"dec10"], start=c(1967,1), frequency=12)
plot(d10, main="CRSP Lower 10% Mothly Returns")

用acf()作时间序列的自相关函数图：

acf()的返回值是一个列表，其中lag相当于k，
acf相当于ρ̂ k。
用plot=FALSE取消默认的图形输出。

ACF图中横轴上下两条水平下是在独立同分布白噪声假设下的加减两倍标准差，
即±2T√。
如果独立同分布白噪声假设成立，
每个ρ̂ k有95%以上的概率落入这两条线之间。

ACF图k=0处总对应ρ̂ 0=1。

上图的ρ̂ 1和ρ̂ 12都超出了界限
（因为是月度数目，横轴的单位是1/12为一个时间点）。
从此图可以认为此投资组合的收益率不是白噪声。

标准库stats的acf()作自相关函数图总是有ρ0=1，
这在其它系数很小时会使得其它系数显得很难分辨，
另外，当时间序列是月度或季度数据时，
横坐标是以年为单位。
所以，forecast包提供了一个类似功能的Acf()函数，如：

TSA包的acf()函数也不做滞后0的点，如：

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast

○○○○○

用单个自相关系数作白噪声检验

如果{Xt}是独立同分布白噪声，
则ρ̂ k(k≥1)近似N(0,1/T)。
若H0是序列为白噪声，
取统计量

t=T‾‾√ρ̂ k

如果|t|>qnorm(1−α/2)，
则拒绝白噪声零假设。
实际中常取α=0.05,
qnorm(1−α/2)≈2，
当ρ̂ 1超出±2/T‾‾√则拒绝H0，
有多个ρ̂ k超出±2/T‾‾√也可拒绝H0，
有一个t统计量值很大（比如超出±3）也可拒绝H0。

在判断{Xt}是否Xt=μ+∑qj=0ψjεt−j
这样的模型时，
根据Bartlet公式，
可取

t=ρ̂ k1T(1+2∑k−1j=1ρ̂ 2j)‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√, k>q

当t超出qnorm(1−α/2)时拒绝这样的模型。

例3.2 这是例3.1 的继续。
有研究者认为小市值股票倾向于在每年的一月份有正的收益率。

为此，用H0:ρ12=0对Ha:ρ12≠0的检验来验证。
如果一月份有取正值的倾向，
则相隔12个月的值会有正相关。

计算ρ̂ 12的值：

tmp1 <- acf(d10, plot=FALSE)
r12 <- tmp1$acf[abs(tmp1$lag-12/12)<1E-10]
r12

## [1] 0.130411

计算t统计量的值，检验p值：

t12 <- sqrt(tmp1$n.used)*r12; t12

## [1] 2.962369

pv <- 2*(1 - pnorm(abs(t12))); pv

## [1] 0.003052812

p值小于0.05,
拒绝了H0:ρ12=0，
这个检验的结果支持一月份效应的存在性。

○○○○○

Ljung-Box白噪声检验

为了检验时间序列样本是否来自白噪声序列，
可以检验ρk=0,k=1,2,…的零假设。
前面检验单个ρk的做法如果针对多个进行检验就有多重检验的第一类错误增大的问题。

Box和Pierce(G. Box and Pierce 1970)提出了混成统计量(Portmanteau statistic)

Q∗(m)=T∑j=1mρ̂ 2j

用来检验

H0:ρ1=⋯=ρm=0⟷Ha:不全为零

在{Xt}是独立白噪声序列条件下，
Q∗(m)近似服从χ2(m)分布。
给定检验水平α，
当Q∗(m)>qchisq(1−α,m)时拒绝H0，
否定白噪声假设。
如果检验的序列是线性时间序列估计的残差序列，
则卡方自由度应改为m减去估计的系数个数。

Ljung和Box(Ljung and Box 1978)对此检验方法进行了改进。
统计量改为

Q(m)=T(T+2)∑j=1mρ̂ 2jT−j

在独立同分布白噪声假设下仍近似服从χ2(m)分布。
当Q(m)>qchisq(1−α,m)时拒绝H0，
否定白噪声假设。
这个检验称为Ljung-Box白噪声检验。
如果检验的序列是线性时间序列估计的残差序列，
则卡方自由度应改为m减去估计的系数个数。
比如，
对ARMA(p,q)模型建模的残差作白噪声检验，
卡方自由度应改为m−(p+q)。

在R软件中，
Box.test(x, type="Ljung-Box")执行Ljung-Box白噪声检验。
Box.test(x, type="Box-Pierce")执行Box-Pierce混成检验。
用fitdf=指定要减去的自由度个数。

R的portes包提供了多种一元和多元白噪声检验功能。

例3.3 检验IBM股票月收益率是否白噪声。

考虑IBM股票从1926-01到2011-09的月度收益率数据，
简单收益率和对数收益率分别考虑。

读入数据：

d <- read_table(
  "m-ibmsp-2611.txt", col_types=cols(
    .default=col_double(),
    date=col_date(format="%Y%m%d")))
ibm <- ts(d[["ibm"]], start=c(1926,1), frequency=12)

读入的是简单收益率的月度数据。
作ACF图：

从ACF来看月度简单收益率是白噪声。

作Ljung-Box白噪声检验，
分别取m=12和m=24：

Box.test(ibm, lag=12, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  ibm
## X-squared = 13.098, df = 12, p-value = 0.362

Box.test(ibm, lag=24, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  ibm
## X-squared = 35.384, df = 24, p-value = 0.0629

在0.05水平下均不拒绝零假设，
支持IBM月度简单收益率是白噪声的零假设。

从简单收益率计算对数收益率，
并进行LB白噪声检验：

Box.test(log(1 + ibm), lag=12, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  log(1 + ibm)
## X-squared = 12.814, df = 12, p-value = 0.3827

Box.test(log(1 + ibm), lag=24, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  log(1 + ibm)
## X-squared = 34.506, df = 24, p-value = 0.07607

在0.05水平下不拒绝零假设。

○○○○○

Box-Pierce检验和Ljung-Box检验受到m取值的影响，
建议采用m≈lnT，
且序列为季度、月度这样的周期序列时，
m应取为周期的整数倍。

例3.4 这是例3.1和例3.2的继续。
对CRSP最低10分位的资产组合的月简单收益率作白噪声检验。

此组合的收益率序列的ACF:

针对m=12和m=24作Ljung-Box白噪声检验：

Box.test(d10, type="Ljung-Box", lag=12)

## 
##  Box-Ljung test
## 
## data:  d10
## X-squared = 41.06, df = 12, p-value = 4.789e-05

Box.test(d10, type="Ljung-Box", lag=24)

## 
##  Box-Ljung test
## 
## data:  d10
## X-squared = 56.246, df = 24, p-value = 0.0002122

在0.05水平下均拒绝零假设，
认为CRSP最低10分位的投资组合的月度简单收益率不是白噪声。

○○○○○

有效市场假设认为收益率是不可预测的，
也就不会有非零的自相关。
但是，股价的决定方式和指数收益率的计算方式等可能会导致在观测到的收益率序列中有自相关性。
高频金融数据中很常见自相关性。

常见的白噪声检验还有TREVOR S. BREUSCH (1978)
和LESLIE G. GODFREY (1978)提出的拉格朗日乘子法检验（LM检验）。
零假设为白噪声，
对立假设为AR、MA或者ARMA。
参见：

• TREVOR S. BREUSCH(1978),
  Testing for Autocorrelation in Dynamic Linear Models,
  Australian Economic Papers 17, pp. 334 – 355
• LESLIE G. GODFREY(1978),
  Testing Against General Autoregressive and Moving Average Error Models When Regressors Include Lagged Dependent Variables,
  Econometrica 46 , S. 1293 – 1302

严平稳

对时间序列{Xt,t∈ℤ}，
如果对任意的t∈ℤ和正整数n, k，
(Xt,…,Xt+n−1)总是与
(Xt+k,…,Xt+n−1+k)同分布，
则称{Xt}为严平稳时间序列。

如果严平稳时间序列{Xt}有二阶矩，
则它也是宽平稳的。

如果宽平稳时间序列{Xt}的任意有限维分布都服从广义正态分布，
则{Xt}也是严平稳的。

如果{εt}为独立同分布零均值白噪声，
方差为σ2，
{ψj}绝对可和或者平方可和，
则线性时间序列

Xt=μ+∑j=0∞ψjεt−j,t∈ℤ

同时宽平稳和严平稳的时间序列。

严平稳遍历性

时间序列建模，
是从理论上无穷长度的随机变量序列{Xt(ω),t=…,−1,0,1,…,ω∈Ω}
中获得一段观测值，
称为一条轨道，
如{X1(ω0),X2(ω0),…,XT(ω0)}，
其中ω0是概率空间Ω中的一个结果点。
虽然T可以越来越大，
但是我们还是仅仅观测到结果空间中的一个点，
希望据此推断所有ω∈Ω的性质。
这需要序列{Xt}具有遍历性。

许多序列没有遍历性。
例如，
Xt≡ξ，
其中ξ为U(0,1)随机变量，
则每条轨道是一个常数，
但不同轨道的常数值不同，
从一条轨道无法推断所有结果。
有如，设随机变量U服从U(0,2π)，
A>0, f∈(0,2π)为常数，则时间序列

Xt=Acos(ft+U), t∈ℤ

是宽平稳时间序列，
但其每条轨道是一条正弦曲线上的离散点，
不同轨道的相位不同，
从一条轨道无法推断所有轨道的分布情况。

如果从时间序列的一条轨道就可以推断出它的所有有限维分布，
就称其为严平稳遍历的。
这里不给出遍历性的严格定义，
仅给出一些严平稳遍历的充分条件。
可以证明，
宽平稳的正态时间序列是严平稳遍历的，
上述由零均值独立同分布白噪声产生的线性序列是严平稳遍历的。

严平稳遍历的性质：

定理:
如果{Xt}是严平稳遍历序列，
则

(1) 强大数律: 如果E|X1|<∞则

limT→∞1T∑t=1TXt=EX1, a.s.

(2) 对任意函数ϕ(x1,x2,…,xm)，
Yt=ϕ(Xt,Xt−1,…,Xt+m−1)也是严平稳遍历序列。