時間序列分析研究的是按時間順序收集的數(shù)據(jù)。相鄰的觀測數(shù)據(jù)通常相互依賴。因此，時間序列分析的技術(shù)需要處理這種相依性。 首先，我們考慮如何在 R 中存儲和處理時間序列。接著，我們處理線性時間序列分析，并展現(xiàn)如何將它用于建模和預(yù)測房屋價格。其次，我們考慮長期趨勢，最后使用協(xié)整的概念來改進(jìn)基本的最小方差對沖比。
數(shù)據(jù)源：http://labfile.oss.aliyuncs.com/courses/882/tarfile.tar.gz
于存儲時間序列數(shù)據(jù)的基本 R 類有 vector、matrix、data.frame 以及 ts 對象。但是，它們可以存儲在這些對象中的數(shù)據(jù)類型相當(dāng)有限。并且，這些表達(dá)方式提供的方法范圍也很有限。不過幸運(yùn)的是，同名的包中的特定對象，zoo、xts 或 timeSeries 對象，對時間序列數(shù)據(jù)提供了更一般的表達(dá)形式。 對每個時間序列分析問題都創(chuàng)建時間序列對象是不必要的，但是復(fù)雜程度較高的分析則需要創(chuàng)建時間序列對象。你可以先將時間序列數(shù)據(jù)存儲成向量形式，再計算數(shù)據(jù)的均值和方差，但如果你想用 decompose 對數(shù)據(jù)做季節(jié)分解，那就必須將數(shù)據(jù)存儲在時間序列對象中。
我們使用蘋果公司股票的日收盤價，創(chuàng)建一個名為 appl 的 zoo 對象，存儲在 CSV 文件 aapl.csv 中。表格的每一行包括一個日期和一個價格，兩項通過逗號分隔。第一行包含了列名（Date 和 Close）。日期格式符合 ISO8601 推薦的基本標(biāo)準(zhǔn)符號（YYYY-MM-DD）。收盤價根據(jù)股票的拆分、股利以及相關(guān)改變進(jìn)行調(diào)整。
> library(zoo)
> aapl <- read.zoo("aapl.csv", sep=",", header=TRUE, format="%Y-%m-%d")
> plot(aapl, main="APPLE Closing Prices on NASDAQ", ylab="Price (USD)", xlab="Date")
使用下面的命令，我們可以提取時間序列開頭部分或結(jié)尾部分。
> head(aapl)
2000-01-03 2000-01-04 2000-01-05 2000-01-06 2000-01-07 2000-01-10
27.58 25.25 25.62 23.40 24.51 24.08
> tail(aapl)
2013-04-17 2013-04-18 2013-04-19 2013-04-22 2013-04-23 2013-04-24
402.80 392.05 390.53 398.67 406.13 405.46
使用下面的命令，可以找出蘋果股價在所有時間中的高點，和這個高點發(fā)生的日期。
> aapl[which.max(aapl)]
2012-09-19
694.86
當(dāng)處理時間序列時，通常收益率更受關(guān)注，價格卻不會。其原因是收益率通常平穩(wěn)。因此我們會計算簡單收益率或連續(xù)復(fù)合收益率（按百分比的形式）。
> ret_simple <- diff(aapl)/lag(aapl, k=-1)*100
> ret_cont <- diff(log(aapl))*100
同時，我們也可以得到簡單收益率的概括統(tǒng)計。在這里，我們使用 coredata 方法來表明我們僅僅關(guān)注股票價格，而非索引（日期）。
> summary(coredata(ret_simple))
Min. 1st Qu. Median Mean 3rd Qu. Max.
-51.85888 -1.32475 0.07901 0.12527 1.55332 13.91131
可以看出，最大的單日損失是 51.86%。我們還可以使用下面的命令獲得這個損失發(fā)生的日期。
> ret_simple[which.min(ret_simple)]
2000-09-29
-51.85888
上網(wǎng)快速搜索可以發(fā)現(xiàn)，這個股價的劇烈變動緣于一個盈利預(yù)警的發(fā)布。我們可以畫出直方圖來加深理解日收益率的相關(guān)頻率。對收益率數(shù)據(jù)進(jìn)行分組時，我們可以使用 break 參數(shù)來指定每組的元素個數(shù)。
> hist(ret_simple, breaks=100, main="Histogram of Simple Returns", xlab="%")
我們也可以把分析限定于時間序列的一個子集（window）中。比如，蘋果股價在 2013 年的最高點可以通過運(yùn)行下面的命令的找到。
> aapl_2013 <- window(aapl, start='2013-01-01', end='2013-12-31')
> aapl_2013[which.max(aapl_2013)]
2013-01-02
545.85
從風(fēng)險管理的角度看，收益率分布的分位數(shù)很有意義。比如，我們使用簡單的歷史模擬法，可以很容易確定一天中置信水平為 99% 的在險價值（Value-at-Risk）。
> quantile(ret_simple, probs=0.01)
1%
-7.042678
因此，在任意給定的一天中，收益率低于 ?7% 的概率只有 1%。但是如果這一天發(fā)生了這樣的情形（每年大約會發(fā)生 2.5 次），7% 將是最小的損失量。 線性時間序列的建模與預(yù)測 線性時間序列的一類重要模型是自回歸單整移動平均（Autoregressive Integrated Moving Average，ARIMA）模型族。ARIMA 模型假定了時間序列的當(dāng)前值只依賴于自身的過去值和某些誤差項的過去值。建立 ARIMA 模型包含了以下 3 個階段。 1．模型識別。 2．模型估計。 3．模型診斷檢驗。 模型識別的階段包括了使用圖方法或信息準(zhǔn)則來確定試驗?zāi)Ｐ偷碾A數(shù)（包含的過去值個數(shù)和過去誤差項個數(shù)）。模型階數(shù)確定之后需要估計模型參數(shù)，通常會使用最小二乘方法或者極大似然方法。最后，為了檢查模型可能存在的缺陷，必須仔細(xì)檢查擬合的模型。這個目的可以通過保證模型殘差的行為符合白噪聲的特點來實現(xiàn)，換句話說，殘差不存在線性依賴。
> hp <- read.zoo("UKHP.csv", sep=",", header=TRUE, format="%Y-%m", FUN=as.yearmon)
參數(shù) FUN 對日期列調(diào)用給定的函數(shù)（as.yearmon，表示月度數(shù)據(jù)點）。為了確認(rèn)指定 as.yearmon 真正存儲了月度數(shù)據(jù)（每個周期 12 個子周期），我們可以查詢數(shù)據(jù)序列的頻率。
> frequency(hp)
[1] 12
個周期（稱為年）中有 12 個子周期（稱為月）。為了深入分析，我們再次計算數(shù)據(jù)的簡單收益率。
> hp_ret <- diff(hp)/lag(hp, k=-1)*100
們使用 forecast 包提供的 auto.arima 函數(shù)，在一步中同時識別最優(yōu)模型并估計參數(shù)。除了收益率序列（hp_ret），函數(shù)還使用了幾個參數(shù)。通過指定 stationary = TRUE，我們將搜索僅僅限于平穩(wěn)模型。同樣地，seasonal = FALSE 將搜索限定于非季節(jié)模式。此外，我們選擇模型時，選擇赤池信息量來度量模型的相對質(zhì)量。
> mod <- auto.arima(hp_ret, stationary=TRUE, seasonal=FALSE, ic="aic")
為了確定擬合系數(shù)的值，我們可以查詢模型的輸出。
> mod
Series: hp_ret
ARIMA(2,0,0) with non-zero mean
Coefficients:
ar1 ar2 mean
0.2299 0.3491 0.4345
s.e. 0.0573 0.0575 0.1519
sigma^2 estimated as 1.118: log likelihood=-390.97
AIC=789.94 AICc=790.1 BIC=804.28
根據(jù)赤池信息量準(zhǔn)則，看起來 AR（2）模型擬合數(shù)據(jù)最好。當(dāng)然，我們也可以使用命令 pacf 畫出偏自相關(guān)函數(shù)，從視覺上來確定滯后階數(shù)。AR（2）模型給出了兩個 AR 系數(shù)、截距（如果模型包含 AR 項，截距實際上是均值）以及相應(yīng)的標(biāo)準(zhǔn)誤。在下面的例子中，因為水平為 5% 的置信區(qū)間沒有包括 0，所以這些統(tǒng)計量都在 5% 的水平上顯著。
> confint(mod)
2.5 % 97.5 %
ar1 0.1174881 0.3422486
ar2 0.2364347 0.4617421
intercept 0.1368785 0.7321623
如果模型包含不顯著的系數(shù)，我們可以使用帶有 fixed 參數(shù)的 arima 函數(shù)重新估計模型，這相當(dāng)于輸入元素為 0 和 NA 的向量。NA 表示相應(yīng)的變量系數(shù)需要估計而 0 表示相應(yīng)的變量系數(shù)需要設(shè)置為 0。
一個快速驗證模型的方法是運(yùn)行下面的命令畫出時間序列的診斷圖。
> tsdiag(mod)
時間序列的診斷 標(biāo)準(zhǔn)化殘差看來沒有表現(xiàn)出波動率聚集，ACF 圖中的殘差自相關(guān)并不顯著，還有自相關(guān)的 Ljung-Box 檢驗的 p 值看起來很高，綜上所以殘差獨立的原假設(shè)不能被拒絕，因此模型看來良好。 為了評估模型對樣本數(shù)據(jù)的擬合良好程度，我們可以畫出原始的月回報（細(xì)的黑色實線）與擬合值（寬的紅色點線）的對比圖形。
> plot(mod$x, lty=1, main="UK house prices: raw data vs. fitted + values", ylab="Return in percent", xlab="Date")
> lines(fitted(mod), lty=2, lwd=2, col="red")
> accuracy(mod)
ME RMSE MAE MPE MAPE MASE
Training set 0.001203096 1.051422 0.8059929 -Inf Inf 0.7154503
ACF1
Training set -0.004848183
這個命令返回平均誤差、均方誤差、平均絕對誤差、平均百分比誤差、平均絕對值百分比誤差和平均絕對比例誤差。
為了預(yù)測接下來 3 個月的月收益率（2013 年 4～6 月），我們使用下面的命令。
> predict(mod, n.ahead=3)
$pred
Apr May Jun
2013 0.5490544 0.7367277 0.5439708
$se
Apr May Jun
2013 1.057401 1.084978 1.165247
所以，我們預(yù)期在接下來的 3 個月中，平均房屋價格稍有增長，但標(biāo)準(zhǔn)誤比較高，大約為 1%。為了畫出帶有標(biāo)準(zhǔn)誤的預(yù)測，我們可以使用下面的命令。
> plot(forecast(mod))
協(xié)整的思想是指尋找非平穩(wěn)時間序列之間的一個線性組合，這個線性組合是一個平穩(wěn)時間序列。因此，協(xié)整方法可以用于檢測非平穩(wěn)時間序列（比如價格）之間的穩(wěn)定長期關(guān)系。 航空燃油的交叉對沖 航空公司很自然需要購買航空燃油。但由于航空燃油價格的波動很劇烈，大部分航空公司會將它們對航空燃油價格變化的風(fēng)險敞口對沖掉一部分。如果市場中缺乏航空燃油 OTC 產(chǎn)品（譯注：OTC 產(chǎn)品指交易所場外柜臺交易產(chǎn)品），航空公司會使用相關(guān)交易所交易的期貨合約（比如，取暖油）來實現(xiàn)對沖。在下面的部分中，我們首先使用經(jīng)典方法導(dǎo)出最優(yōu)對沖比率，這種方法僅僅考慮兩種價格之間短期波動。然后考察價格之間的長期穩(wěn)定聯(lián)系，進(jìn)而改進(jìn)最優(yōu)對沖比。
> library(urca)
Warning message:
程輯包‘urca’是用R版本3.6.3 來建造的
> prices <- read.zoo("JetFuelHedging.csv", sep=",", FUN=as.yearmon, format="%Y-%m", header=TRUE)
現(xiàn)在我們僅僅考慮兩種商品的短期行為（月價格改變）。通過擬合一個用取暖油價格變化來解釋航空燃油價格的線性模型，可以推導(dǎo)出兩種商品的最小方差對沖比率。線性模型中的 beta 系數(shù)就是最優(yōu)對沖比。
> simple_mod <- lm(diff(prices$JetFuel) ~ diff(prices$HeatingOil) + 0)
函數(shù) lm（用于線性模型）估計了航空燃油價格變動對取暖油價格變動的最佳擬合系數(shù)。+0 項表示截距設(shè)置為 0，意味著航空公司不持有現(xiàn)金。
> summary(simple_mod)
Call:
lm(formula = diff(prices$JetFuel) ~ diff(prices$HeatingOil) +
0)
Residuals:
Min 1Q Median 3Q Max
-0.52503 -0.02968 0.00131 0.03237 0.39602
Coefficients:
Estimate Std. Error t value Pr(>|t|)
diff(prices$HeatingOil) 0.89059 0.03983 22.36 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0846 on 189 degrees of freedom
Multiple R-squared: 0.7257, Adjusted R-squared: 0.7242
F-statistic: 499.9 on 1 and 189 DF, p-value: < 2.2e-16
結(jié)果得到了對沖比率為 0.89059 和殘差標(biāo)準(zhǔn)差為 0.0846。但是，交叉對沖并非完美無缺，推導(dǎo)出的對沖組合結(jié)果仍然有風(fēng)險。 現(xiàn)在，我們通過考察航空燃油和取暖油期貨價格之間存在的長期關(guān)系，嘗試改進(jìn)方差比。
> plot(prices$JetFuel, main="Jet Fuel and Heating Oil Prices", xlab="Date", ylab="USD")
> lines(prices$HeatingOil, col="red")
我們使用 Engle 和 Granger 的兩步估計方法。首先，使用增強(qiáng)的 Dickey-Fuller 檢驗（augmented Dickey-Fuller test，ADF 檢驗）對兩個序列進(jìn)行單位根檢驗（非平穩(wěn)性）。
> jf_adf <- ur.df(prices$JetFuel, type="drift")
> summary(jf_adf)
###############################################
# 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
-1.06212 -0.05015 0.00566 0.07922 0.38086
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.03050 0.02177 1.401 0.16283
z.lag.1 -0.01441 0.01271 -1.134 0.25845
z.diff.lag 0.19471 0.07250 2.686 0.00789 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.159 on 186 degrees of freedom
Multiple R-squared: 0.04099, Adjusted R-squared: 0.03067
F-statistic: 3.975 on 2 and 186 DF, p-value: 0.0204
Value of test-statistic is: -1.1335 0.9865
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
結(jié)果顯示，因為檢驗統(tǒng)計量值 ?1.1335 大于臨界值 ?3.46，所以在 1% 的置信水平上不能拒絕非平穩(wěn)（航空燃油時間序列包含一個單位根）的原假設(shè)。同樣的結(jié)果對取暖油也成立（檢驗統(tǒng)計量是 ?1.041）。
> ho_adf <- ur.df(prices$HeatingOil, type="drift")
> summary(ho_adf)
###############################################
# 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.78839 -0.06344 -0.00128 0.07418 0.49186
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.02815 0.02115 1.331 0.1847
z.lag.1 -0.01314 0.01263 -1.041 0.2992
z.diff.lag 0.14419 0.07296 1.976 0.0496 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1534 on 186 degrees of freedom
Multiple R-squared: 0.02418, Adjusted R-squared: 0.01369
F-statistic: 2.304 on 2 and 186 DF, p-value: 0.1027
Value of test-statistic is: -1.041 0.9002
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
現(xiàn)在我們可以繼續(xù)估計靜態(tài)均衡模型，并使用 ADF 方法檢驗時間序列的殘差是否平穩(wěn)。請注意，目前的研究序列是上一步的估計結(jié)果，因此，我們現(xiàn)在必須使用不同的臨界值。
> mod_static <- summary(lm(prices$JetFuel ~ prices$HeatingOil))
> error <- residuals(mod_static)
> error_cadf <- ur.df(error, type="none")
> summary(error_cadf)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.19162 -0.03385 -0.00516 0.02879 0.47426
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.74476 0.08357 -8.912 4.46e-16 ***
z.diff.lag 0.12304 0.07264 1.694 0.092 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07107 on 187 degrees of freedom
Multiple R-squared: 0.3417, Adjusted R-squared: 0.3347
F-statistic: 48.53 on 2 and 187 DF, p-value: < 2.2e-16
Value of test-statistic is: -8.912
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
得到的檢驗統(tǒng)計量是 ?8.912，而規(guī)模為 200 的樣本在 1% 的置信水平上的臨界值為 ?2.85。所以，我們拒絕非平穩(wěn)的原假設(shè)。因此我們發(fā)現(xiàn)了兩個協(xié)整變量并且可以進(jìn)行第二步，兩個協(xié)整變量意味著一個誤差修正模型（ECM）。ECM 是一個動態(tài)模型，刻畫了系統(tǒng)如何（以及多快）返回之前估計的靜態(tài)均衡，這個靜態(tài)均衡存儲在變量 mod_static 中。
> djf <- diff(prices$JetFuel)
> dho <- diff(prices$HeatingOil)
> error_lag <- lag(error, k=-1)
> mod_ecm <- lm(djf ~ dho + error_lag)
> summary(mod_ecm)
Call:
lm(formula = djf ~ dho + error_lag)
Residuals:
Min 1Q Median 3Q Max
-0.19332 -0.03391 -0.00101 0.02128 0.44942
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.001511 0.005013 0.301 0.763
dho 0.899489 0.032547 27.637 <2e-16 ***
error_lag -0.655301 0.066303 -9.883 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06892 on 187 degrees of freedom
Multiple R-squared: 0.8189, Adjusted R-squared: 0.817
F-statistic: 422.9 on 2 and 187 DF, p-value: < 2.2e-16
通過考察航空燃油價格和取暖油價格之間存在的長期聯(lián)系（協(xié)整），對沖比現(xiàn)在稍微提高了一點（0.89948），并且殘差標(biāo)準(zhǔn)差顯著地降低了（?0.65530）。誤差項的系數(shù)為負(fù)（?0.65530）：兩個價格之間原先較大的偏差有所修正，價格向它們的長期穩(wěn)定關(guān)系移動。
