【摘要】时间序列数据是指为单个实体随时间变化而收集的数据。
◾Cross-section data
◾Panel data
使用时间序列数据的两个应用
◾Dynamic casusl effect: tax on beer on future fatality rates
◾Forcasting: temperature, GDP, unemployment rate

时间序列（time series）是随机变量Y1、Y2、……Yt的一个序列，它是由等距的时间点序列索引的。
一个时间序列的均值函数就是该时间序列在某个时间索引t上的期望值。一般情况下，某个时间序列在某个时间索引t1的均值并不等于该时间序列在另一个不同的时间索引t2的均值。
自协方差函数及自相关函数是衡量构成时间序列的随机变量在不同时间点上相互线性依赖性的两个重要函数。自相关函数通常缩略为ACF函数。ACF函数是对称的，但是无单位，其绝对值被数值1约束，即当两个时间序列索引之间的自相关度是1或-1,就代表两者之间存在完全线性依赖或相关，而当相关度是0时，就代表完全线性无关。
平稳性：实质描述的是一个时间序列的概率表现不会随着时间的流逝而改变。常用的平稳性的性质有严格平稳和弱平稳两个版本。tseries包的adf.test()函数可以检验时间序列的平稳性，返回的p值小于0.05则表示是平稳的。
白噪声是一个平稳过程，因为它的均值和方差都是常数。
随机漫步的均值是常数（不带漂移的随机漫步），但它的方差是随着时间的变化而不同的，因此它是不平稳的。
自回归模型（Autoregressive models， AR）来源于要让一个简单模型根据过去有限窗口时间里的最近值来解释某个时间序列当前值的想法。
自回归条件异方差模型：ARIMA模型的关键前提条件是，虽然序列本身是非平稳的，但是我们可以运用某个变换来获得一个平稳的序列。像这样为非平稳时间序列构建模型的方法之一是作出一个假设，假设该模型非平稳的原因是该模型的方差会以一种可预见的方式随时间变化，这样就可以把方差随时间的变化建模为一个自回归过程，这种模型被称为自回归条件异方差模型（ARCH）。加入了移动平均方差成分的ARCH模型称为广义自回归条件异方差模型（GARCH）。

zoo是一个R语言类库，zoo类库中定义了一个名为zoo的S3类型对象，用于描述规则的和不规则的有序的时间序列数据。zoo对象是一个独立的对象，包括索引、日期、时间，只依赖于基础的R环境，zooreg对象继承了zoo对象，只能用于规则的的时间序列数据。
◾zoo对象包括两部分组成，数据部分、索引部分。R语言时间序列基础库zoo
函数定义：
zoo(x = NULL, order.by = index(x), frequency = NULL)
参数列表：
x: 数据部分，允许向量，矩阵，因子
order.by: 索引部分，唯一字段，用于排序
frequency: 每个时间单元显示的数量
数据显示格式
以年+季度格式输出 as.yearqtr()
以年+月份格式输出as.yearmon()
◾xts是对时间序列数据(zoo)的一种扩展实现，目标是为了统一时间序列的操作接口。可扩展的时间序列xts

pacman::p_load(AER,dynlm,forecast,readxl,stargazer,scales,quantmod,urca)

◾时间序列数据和序列相关性

# 包'quantmod'对时间序列数据进行画图和计算
library(quantmod)
USMacroSWQ <- read_xlsx("C:/Users/apple/Desktop/us_macro_quarterly.xlsx",
                        sheet = 1,
                        col_types = c("text", rep("numeric", 9)))
#col_types = c("text", rep("numeric", 9))第一列是文本型，其余都是数字型

# 日期列格式化，yearqtr()以年季度显示时间
USMacroSWQ$...1 <- as.yearqtr(USMacroSWQ$...1, format = "%Y:0%q")
# 调整列名
colnames(USMacroSWQ) <- c("Date", "GDPC96", "JAPAN_IP", "PCECTPI",
                          "GS10", "GS1", "TB3MS", "UNRATE", "EXUSUK", "CPIAUCSL")

# 把GDP定义为时间序列变量,xts
GDP <- xts(USMacroSWQ$GDPC96, USMacroSWQ$Date)["1960::2013"]
# 把GDP growth定义为时间序列变量，近似计算log(GDP/lag(GDP))~x2/x1-1
GDPGrowth <- xts(400 * log(GDP/lag(GDP)))

#as.zoo：GDP定义为时间序列变量，不用定义x/y
plot(log(as.zoo(GDP)),
     col = "steelblue",
     lwd = 2,
     ylab = "Logarithm",
     xlab = "Date",
     main = "U.S. Quarterly Real GDP")

plot(as.zoo(GDPGrowth),
     col = "steelblue",
     lwd = 2,
     ylab = "Logarithm",
     xlab = "Date",
     main = "U.S. Real GDP Growth Rates")

符号、标签、差异、对数法和增长率

#lag(x, k = 1, ...)k=1选取前一期数据
#diff(Y)差分，增长率
#计算对数、年增长率和增长率
quants <- function(series) {
  s <- series
  return(
    data.frame("Level" = s,
               "Logarithm" = log(s),
               "AnnualGrowthRate" = 400 * log(s / lag(s)),
               "1stLagAnnualGrowthRate" = lag(400 * log(s / lag(s))))
  ) }
quants(GDP["2011-07::2013-01"])

> quants(GDP["2011-07::2013-01"])
           Level Logarithm AnnualGrowthRate X1stLagAnnualGrowthRate
2011 Q3 15062.14  9.619940               NA                      NA
2011 Q4 15242.14  9.631819        4.7518062                      NA
2012 Q1 15381.56  9.640925        3.6422231               4.7518062
2012 Q2 15427.67  9.643918        1.1972004               3.6422231
2012 Q3 15533.99  9.650785        2.7470216               1.1972004
2012 Q4 15539.63  9.651149        0.1452808               2.7470216
2013 Q1 15583.95  9.653997        1.1392015               0.1452808

自相关(Autocorrelation)

#计算GDPGrowth的前四个样本自相关性
acf(na.omit(GDPGrowth), lag.max = 4, plot = F)
#绘制自相关图acf(x,lag=):lag表示延迟阶数.若用户不特殊指定的话，系统会根据序列长度自动指定延迟阶数.
>Autocorrelations of series ‘na.omit(GDPGrowth)’, by lag

 0.00  0.25  0.50  0.75  1.00 
1.000 0.352 0.273 0.114 0.106 

经济时间序列的其他例子

# define series as xts objects
USUnemp <- xts(USMacroSWQ$UNRATE, USMacroSWQ$Date)["1960::2013"]
DollarPoundFX <- xts(USMacroSWQ$EXUSUK, USMacroSWQ$Date)["1960::2013"]
JPIndProd <- xts(log(USMacroSWQ$JAPAN_IP), USMacroSWQ$Date)["1960::2013"]
# attach NYSESW data
data("NYSESW")
NYSESW <- xts(Delt(NYSESW))
# divide plotting area into 2x2 matrix
par(mfrow = c(2, 2))
# plot the series
plot(as.zoo(USUnemp),
     col = "steelblue",
     lwd = 2,
     ylab = "Percent",
     xlab = "Date",
     main = "US Unemployment Rate",
     cex.main = 1)
plot(as.zoo(DollarPoundFX),
     col = "steelblue",
     lwd = 2,
     ylab = "Dollar per pound",
     xlab = "Date",
     main = "U.S. Dollar / B. Pound Exchange Rate",
     cex.main = 1)
plot(as.zoo(JPIndProd),
     col = "steelblue",
     lwd = 2,
     ylab = "Logarithm",
     xlab = "Date",
     main = "Japanese Industrial Production",
     cex.main = 1)
plot(as.zoo(NYSESW),
     col = "steelblue",
     lwd = 2,
     ylab = "Percent per Day",
     xlab = "Date",
     main = "New York Stock Exchange Composite Index",
     cex.main = 1)

◾失业率在经济衰退期间上升，在经济复苏和增长期间下降。
◾美元/磅的汇率显示了一个确定性的模式，直到布雷顿森林系统的结束。
◾日本的工业生产呈上升趋势和下降趋势。
◾在纽约证券交易所综合指数的每日变化似乎在零线附近随机波动。

◾自回归

自回归模型被大量用于经济预测。自回归模型将时间序列变量与其过去的值联系起来。

一阶自回归模型(The First-Order Autoregressive Model)

#数据子集
GDPGRSub <- GDPGrowth["1962::2012"]
# 建立模型
# ar.ols自回归模型,与lm()一样
# order.max = 1一阶自回归
ar.ols(GDPGRSub, 
       order.max = 1, 
       demean = F, 
       intercept = T)

#用lm()建立模型
# length of data set
N <-length(GDPGRSub)
GDPGR_level <- as.numeric(GDPGRSub[-1])
GDPGR_lags <- as.numeric(GDPGRSub[-N])
# estimate the model
armod <- lm(GDPGR_level ~ GDPGR_lags)
armod
# 对估计回归系数的稳健summary
coeftest(armod, vcov. = vcovHC, type = "HC1")

> coeftest(armod, vcov. = vcovHC, type = "HC1")

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 1.994986   0.351274  5.6793 4.691e-08 ***
GDPGR_lags  0.338436   0.076188  4.4421 1.470e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

对GDP增长的应用

library(forecast)
# assign GDP growth rate in 2012:Q4
new <- data.frame("GDPGR_lags" = GDPGR_level[N-1])
# forecast GDP growth rate in 2013:Q1
forecast(armod, newdata = new)

> forecast(armod, newdata = new)
  Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
1       2.044155 -2.036225 6.124534 -4.213414 8.301723

p阶的自回归模型

# 建立AR（2）估计模型,动态滞后模型，加载 dynlm 包。
GDPGR_AR2 <- dynlm(ts(GDPGR_level) ~ L(ts(GDPGR_level)) + L(ts(GDPGR_level), 2))
coeftest(GDPGR_AR2, vcov. = sandwich)
# Rˆ2
summary(GDPGR_AR2)$r.squared
# SER
summary(GDPGR_AR2)$sigma
# AR(2) 预测2013:Q1GDP growth 
forecast <- c("2013:Q1" = coef(GDPGR_AR2) %*% c(1, GDPGR_level[N-1], GDPGR_level[N-2]))

> coeftest(GDPGR_AR2, vcov. = sandwich)

t test of coefficients:

                      Estimate Std. Error t value  Pr(>|t|)    
(Intercept)           1.631747   0.402023  4.0588 7.096e-05 ***
L(ts(GDPGR_level))    0.277787   0.079250  3.5052 0.0005643 ***
L(ts(GDPGR_level), 2) 0.179269   0.079951  2.2422 0.0260560 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


> summary(GDPGR_AR2)$r.squared
[1] 0.1425484


> summary(GDPGR_AR2)$sigma
[1] 3.132122

能打败市场吗？

◾有效资本市场的理论表明，股票价格体现了目前所有可用的信息。
◾如果这一假设成立，就不可能利用过去收益的公开信息建立一个有用模型去估计预测未来股票收益：如果有可能预测市场，交易员将能够套利，例如，通过依赖AR（(2)）模型，他们将使用尚未定价的信息来推动价格，直到预期收益为零。

SReturns <- read_xlsx("C:/Users/apple/Desktop/Stock_Returns_1931_2002.xlsx",
                      sheet = 1,
                      col_types = "numeric")
# 把第三、四列定义为时间序列数据
StockReturns <- ts(SReturns[, 3:4],
                   start = c(1931, 1),
                   end = c(2002, 12),
                   frequency = 12)

# 建立AR模型:
# AR(1)
SR_AR1 <- dynlm(ExReturn ~ L(ExReturn),
                data = StockReturns, start = c(1960, 1), end = c(2002, 12))
#summary(SR_AR1)
# AR(2)
SR_AR2 <- dynlm(ExReturn ~ L(ExReturn) + L(ExReturn, 2),
                data = StockReturns, start = c(1960, 1), end = c(2002, 12))
#summary(SR_AR2)
# AR(4)
SR_AR4 <- dynlm(ExReturn ~ L(ExReturn) + L(ExReturn, 1:4),
                data = StockReturns, start = c(1960, 1), end = c(2002, 12))
#summary(SR_AR4)
# 计算稳健标准误
rob_se <- list(sqrt(diag(sandwich(SR_AR1))),
               sqrt(diag(sandwich(SR_AR2))),
               sqrt(diag(sandwich(SR_AR4))))

# 用stargazer()生成表格
stargazer(SR_AR1, SR_AR2, SR_AR4,
          title = "Autoregressive Models of Monthly Excess Stock Returns",
          header = FALSE,
          model.numbers = F,
          omit.table.layout = "n",
          digits = 3,
          column.labels = c("AR(1)", "AR(2)", "AR(4)"),
          dep.var.caption = "Dependent Variable: Excess Returns on the CSRP Value-Weighted Index",
          dep.var.labels.include = FALSE,
          covariate.labels = c("$excess return_{t-1}$", "$excess return_{t-2}$",
                               "$excess return_{t-3}$", "$excess return_{t-4}$",
                               "Intercept"),
          se = rob_se,
          omit.stat = "rsq")

◾研究结果与有效金融市场的假设一致：在任何估计模型中都没有统计上显著的系数，所有系数都为零的假设不能被拒绝。
◾barR^2在所有型号中几乎为零，在AR(4)模型中几乎为负值。这表明这些模型都没有对预测股票回报有用。

其他预测和ADL模型(autoregressive distributed lag )

利用期限利差预测GDP增长

◾长期和短期国债的利率与宏观经济状况密切相关。虽然这两种债券的利率都有相同的长期趋势，但在短期内它们的表现完全不同。
◾两种到期日不同的债券的利率差值术语称为期限利差。

# 3月期国债利率
TB3MS <- xts(USMacroSWQ$TB3MS, USMacroSWQ$Date)["1960::2012"]
# 10年期国债利率
TB10YS <- xts(USMacroSWQ$GS10, USMacroSWQ$Date)["1960::2012"]
# 期限差
TSpread <- TB10YS - TB3MS
# reproduce Figure 14.2 (a) of the book
plot(merge(as.zoo(TB3MS), as.zoo(TB10YS)),
     plot.type = "single",
     col = c("darkred", "steelblue"),
     lwd = 2,
     xlab = "Date",
     ylab = "Percent per annum",
     main = "Interest Rates")

#定义将年份转换为yearqtr类的函数
YToYQTR <- function(years) {
  return(
    sort(as.yearqtr(sapply(years, paste, c("Q1", "Q2", "Q3", "Q4"))))
  ) }

# 经济衰退
recessions <- YToYQTR(c(1961:1962, 1970, 1974:1975, 1980:1982, 1990:1991, 2001, 2007:2008))

# 为经济衰退添加颜色阴影
xblocks(time(as.zoo(TB3MS)),
        c(time(TB3MS) %in% recessions),
        col = alpha("steelblue", alpha = 0.3))
# 添加图例
legend("topright",
       legend = c("TB3MS", "TB10YS"),
       col = c("darkred", "steelblue"),
       lwd = c(2, 2))

# reproduce Figure 14.2 (b) of the book
plot(as.zoo(TSpread),
     col = "steelblue",
     lwd = 2,
     xlab = "Date",
     ylab = "Percent per annum",
     main = "Term Spread")
# 为经济衰退添加颜色阴影
xblocks(time(as.zoo(TB3MS)),
        c(time(TB3MS) %in% recessions),
        col = alpha("steelblue", alpha = 0.3))

◾在衰退前，长期债券和短期票据之间的差距缩小，因此利差急剧下降至零，在经济压力时期甚至变为负。
◾此信息可用于改善对未来的GDP增长预测。

# 转换为时间序列数据
GDPGrowth_ts <- ts(GDPGrowth,
                   start = c(1960, 1),
                   end = c(2013, 4),
                   frequency = 4)
TSpread_ts <- ts(TSpread,
                 start = c(1960, 1),
                 end = c(2012, 4),
                 frequency = 4)
# 联合两个ts对象
ADLdata <- ts.union(GDPGrowth_ts, TSpread_ts)

# 建立GDP growth的ADL(2,1)模型 
GDPGR_ADL21 <- dynlm(GDPGrowth_ts ~ L(GDPGrowth_ts) + L(GDPGrowth_ts, 2) + L(TSpread_ts),
                     start = c(1962, 1), end = c(2012, 4))
coeftest(GDPGR_ADL21, vcov. = sandwich)
# 2012:Q3 / 2012:Q4 GDP growth 和 term spread数据
subset <- window(ADLdata, c(2012, 3), c(2012, 4))
# ADL(2,1)模型预测 2013:Q1GDP growth 
ADL21_forecast <- coef(GDPGR_ADL21) %*% c(1, subset[2, 1], subset[1, 1], subset[2, 2])
ADL21_forecast
# 计算预测错误
window(GDPGrowth_ts, c(2013, 1), c(2013, 1)) - ADL21_forecast

# 建立GDP growth的ADL(2,2)模型 
GDPGR_ADL22 <- dynlm(GDPGrowth_ts ~ L(GDPGrowth_ts) + L(GDPGrowth_ts, 2) 
                     + L(TSpread_ts) + L(TSpread_ts, 2), 
                     start = c(1962, 1), end = c(2012, 4))
coeftest(GDPGR_ADL22, vcov. = sandwich)
# ADL(2,2)模型预测 2013:Q1GDP growth
ADL22_forecast <- coef(GDPGR_ADL22) %*% c(1, subset[2, 1], subset[1, 1], subset[2, 2], subset[1, 2])
ADL22_forecast
# 计算预测错误
window(GDPGrowth_ts, c(2013, 1), c(2013, 1)) - ADL22_forecast

# 比较adj. R2
c("Adj.R2 AR(2)" = summary(GDPGR_AR2)$r.squared,
  "Adj.R2 ADL(2,1)" = summary(GDPGR_ADL21)$r.squared,
  "Adj.R2 ADL(2,2)" = summary(GDPGR_ADL22)$r.squared)

# term spread系数的F检验
linearHypothesis(GDPGR_ADL22, 
                 c("L(TSpread_ts)=0", "L(TSpread_ts, 2)=0"),
                 vcov. = sandwich)

> ADL21_forecast
         [,1]
[1,] 2.241689


> window(GDPGrowth_ts, c(2013, 1), c(2013, 1)) - ADL21_forecast
          Qtr1
2013 -1.102487


> ADL22_forecast
         [,1]
[1,] 2.274407


> window(GDPGrowth_ts, c(2013, 1), c(2013, 1)) - ADL22_forecast
          Qtr1
2013 -1.135206

>    Adj.R2 AR(2) Adj.R2 ADL(2,1) Adj.R2 ADL(2,2) 
      0.1425484       0.1743996       0.1855245 


>Linear hypothesis test

Hypothesis:
L(TSpread_ts) = 0
L(TSpread_ts, 2) = 0

Model 1: restricted model
Model 2: GDPGrowth_ts ~ L(GDPGrowth_ts) + L(GDPGrowth_ts, 2) + L(TSpread_ts) + 
    L(TSpread_ts, 2)

Note: Coefficient covariance matrix supplied.

  Res.Df Df      F  Pr(>F)  
1    201                    
2    199  2 4.4344 0.01306 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

例子：一个趋势的随机行走模型 random walks

arima.sim函数可以拟合平稳AR序列、MA序列、平稳ARMA序列以及ARIMA序列。
arima.sim(n, list(ar=,ma=,order=),sd=)
n:拟合序列长度
list：指定具体参数模型参数，其中：
（1）拟合平稳AR（p）模型，要给出自回归系数，如果指定拟合AR模型为非平稳模型，系统会报错。
（2）拟合MA(q)模型，要给出移动平均系数。
（3）拟合平稳ARMA(p,q)模型，需要同时给出自回归系数和移动平均系数，如果指定模型为非平稳系统会报错。
（4）拟合ARIMA（p，d，q）模型，除了需要给出自回归系数和移动平均系数，还需要增加order选项，order=c（p，d，q），其中p为自回归系数，d为差分阶数，q为移动平均数。
sd：指定序列的标准差，默认sd=1.

# simulate and plot random walks starting at 0
set.seed(1)
RWs <- ts(replicate(n = 4, 
                    arima.sim(model = list(order = c(0, 1 ,0)), n = 100)))

matplot(RWs, 
        type = "l", 
        col = c("steelblue", "darkgreen", "darkred", "orange"), 
        lty = 1, 
        lwd = 2,
        main = "Four Random Walks",
        xlab = "Time",
        ylab = "Value")

# simulate and plot random walks with drift偏移
set.seed(1)

RWsd <- ts(replicate(n = 4, 
                     arima.sim(model = list(order = c(0, 1, 0)), 
                               n = 100,
                               mean = -0.2)))

matplot(RWsd, 
        type = "l", 
        col = c("steelblue", "darkgreen", "darkred", "orange"), 
        lty = 1, 
        lwd = 2,
        main = "Four Random Walks with Drift",
        xlab = "Time",
        ylab = "Value")

具有偏移的随机行走模型，允许模拟一系列向上或向下移动的趋势。如果beta_0为正，则该系列会向上偏移，如果beta_0为负，则遵循向下趋势。

random walks starting at 0

random walks with drift偏移