R語言COPULA和金融時間序列案例

原文 http://tecdat.cn/?p=3385

最近我被要求撰寫關於金融時間序列的copulas的調查。 從讀取數據中得到各類模型的描述，包括一些圖形和統計輸出。

> temp < - tempfile（）
> download.file（
+「http://freakonometrics.free.fr/oil.xls",temp）

> oil = read.xlsx（temp，sheetName =「DATA」，dec =「，」）
> oil = read.xlsx（「D：\\ home \\ acharpen \\ mes documents \\ oil.xls」，sheetName =「DATA」）

而後咱們能夠繪製這三個時間序列

1 1997-01-10 2.73672 2.25465 3.3673 1.5400

2 1997-01-17 -3.40326 -6.01433 -3.8249 -4.1076

3 1997-01-24 -4.09531 -1.43076 -6.6375 -4.6166

4 1997-01-31 -0.65789 0.34873 0.7326 -1.5122

5 1997-02-07 -3.14293 -1.97765 -0.7326 -1.8798

6 1997-02-14 -5.60321 -7.84534 -7.6372 -11.0549

​

這個想法是在這裏使用一些多變量ARMA-GARCH過程。這裏的啓發式是第一部分用於模擬時間序列平均值的動態，第二部分用於模擬時間序列方差的動態。本文考慮了兩種模型

• 關於ARMA模型殘差的多變量GARCH過程（或方差矩陣動力學模型）
• 關於ARMA-GARCH過程殘差的多變量模型（基於copula）

所以，這裏將考慮不一樣的系列，做爲不一樣模型的殘差得到。咱們還能夠將這些殘差標準化。來自ARMA模型

> fit1 = arima（x = dat [，1]，order = c（2,0,1））
> fit2 = arima（x = dat [，2]，order = c（1,0,1））
> fit3 = arima（x = dat [，3]，order = c（1,0,1））
> m < - apply（dat_arma，2，mean）
> v < - apply（dat_arma，2，var）
> dat_arma_std < - t（（t（dat_arma）-m）/ sqrt（v））

或者來自ARMA-GARCH模型

> fit1 = garchFit（formula = ~arma（2,1）+ garch（1,1），data = dat [，1]，cond.dist =「std」）
> fit2 = garchFit（formula = ~arma（1,1）+ garch（1,1），data = dat [，2]，cond.dist =「std」）
> fit3 = garchFit（formula = ~arma（1,1）+ garch（1,1），data = dat [，3]，cond.dist =「std」）
> m_res < - apply（dat_res，2，mean）
> v_res < - apply（dat_res，2，var）
> dat_res_std = cbind（（dat_res [，1] -m_res [1]）/ sqrt（v_res [1]），（dat_res [，2] -m_res [2]）/ sqrt（v_res [2]），（dat_res [ ，3] -m_res [3]）/ SQRT（v_res [3]））

===

多變量GARCH模型

能夠考慮的第一個模型是協方差矩陣的多變量EWMA，

> ewma = EWMAvol（dat_res_std，lambda = 0.96）

要想象波動性，請使用

> emwa_series_vol = function（i = 1）{
+ lines（Time，dat_arma [，i] + 40，col =「gray」）
+ j = 1
+ if（i == 2）j = 5
+ if（i == 3）j = 9

​

隱含的相關性在這裏

> emwa_series_cor = function（i = 1，j = 2）{
+ if（（min（i，j）== 1）＆（max（i，j）== 2））{
+ a = 1; B = 9; AB = 3}
+ r = ewma $ Sigma.t [，ab] / sqrt（ewma $ Sigma.t [，a] *
+ ewma $ Sigma.t [，b]）
+ plot（Time，r，type =「l」，ylim = c（0,1））
+}

​

也可能有一些多變量GARCH，即BEKK（1,1）模型，例如使用：

> bekk = BEKK11（dat_arma）
> bekk_series_vol function（i = 1）{
+ plot（Time， $ Sigma.t [，1]，type =「l」，
+ ylab = （dat）[i]，col =「white」，ylim = c（0,80））
+ lines（Time，dat_arma [，i] + 40，col =「gray」）
+ j = 1
+ if（i == 2）j = 5

+ if（i == 3）j = 9

> bekk_series_cor = function（i = 1，j = 2）{
+ a = 1; B = 5; AB = 2}
+ a = 1; B = 9; AB = 3}
+ a = 5; B = 9; AB = 6}
+ r = bk $ Sigma.t [，ab] / sqrt（bk $ Sigma.t [，a] *
+ bk $ Sigma.t [，b]）

​

從單變量GARCH模型中模擬殘差

第一步多是考慮殘差的一些靜態（聯合）分佈。單變量邊際分佈是

​

而關節密度的輪廓（使用雙變量核估計器得到） 

​

也能夠將copula密度可視化（頂部有一些非參數估計，下面是參數copula）

> copula_NP = function（i = 1，j = 2）{
+ n = nrow（uv）
+ s = 0.3

+ norm.cop < - normalCopula（0.5）
+ norm.cop < - normalCopula（fitCopula（norm.cop，uv）@estimate）
+ dc = function（x，y）dCopula（cbind（x，y），norm.cop）


+ ylab = names（dat）[j]，zlab =「copule Gaussienne」，ticktype =「detailed」，zlim = zl）
+
+ t.cop < - tCopula（0.5，df = 3）
+ t.cop < - tCopula（t.fit [1]，df = t.fit [2]）


+ ylab = names（dat）[j]，zlab =「copule de Student」，ticktype =「detailed」，zlim = zl）
+}

​

能夠考慮這個​功能，

​

計算三對系列的經驗版本，並將其與一些參數版本進行比較，

>

> lambda = function（C）{
+ l = function（u）pcopula（C，cbind（u，u））/ u
+ v = Vectorize（l）（u）
+ return（c（v，rev（v）））
+}
>

> graph_lambda = function（i，j）{
+ X = dat_res
+ U = rank（X [，i]）/（nrow（X）+1）
+ V = rank（X [，j]）/（nrow（X）+1）

+ normal.cop < - normalCopula（.5，dim = 2）
+ t.cop < - tCopula（.5，dim = 2，df = 3）
+ fit1 = fitCopula（normal.cop，cbind（U，V），method =「ml」）
d（U，V），method =「ml」）
+ C1 = normalCopula（fit1 @ copula @ parameters，dim = 2）
+ C2 = tCopula（fit2 @ copula @ parameters [1]，dim = 2，df = trunc（fit2 @ copula @ parameters [2]））
+

​

但人們可能想知道相關性是否隨時間穩定。E：

> time_varying_correl_2 = function（i = 1，j = 2，
+ nom_arg =「Pearson」）{
+ uv = dat_arma [，c（i，j）]
nom_arg））[1,2]
+}
> time_varying_correl_2（1,2）

> time_varying_correl_2（1,2，「spearman」）

> time_varying_correl_2（1,2，「kendall」）

​

斯皮爾曼與時間排名相關

​

或肯德爾 

​

爲了模型的相關性，考慮DCC模型（S）

> m2 = dccFit（dat_res_std）
> m3 = dccFit（dat_res_std，type =「Engle」）
> R2 = m2 $ rho.t
> R3 = m3 $ rho.t

​

要得到一些預測， 使用例如

> garch11.spec = ugarchspec（mean.model = list（armaOrder = c（2,1）），variance.model = list（garchOrder = c（1,1），model =「GARCH」））
> dcc.garch11.spec = dccspec（uspec = multispec（replicate（3，garch11.spec）），dccOrder = c（1,1），
distribution =「mvnorm」）
> dcc.fit = dccfit（dcc.garch11.spec，data = dat）
> fcst = dccforecast（dcc.fit，n.ahead = 200）

​