#### Selected R commands and outputs for Chapter 7 #### Output is edited to remove irrelevant parts. #### > require(MTS) > zt=matrix(rnorm(2000),400,5) ### Another set of random noises > MarchTest(zt) #### Multivariate volatility tests Q(m) of squared series(LM test): Test statistic: 5.054594 p-value: 0.8875013 Rank-based Test: Test statistic: 9.377031 p-value: 0.4967344 Q_k(m) of squared series: Test statistic: 223.293 p-value: 0.8867984 Robust Test(5%) : 269.4988 p-value: 0.1894433 ### IBM and SP log returns > da=read.table("m-ibmsp-6111.txt",header=T) > rtn=log(da[,2:3]+1) > at=scale(rtn,scale=F) ## Remove sample means > MarchTest(at) Q(m) of squared series(LM test): Test statistic: 38.06663 p-value: 3.695138e-05 Rank-based Test: Test statistic: 108.3798 p-value: 0 Q_k(m) of squared series: Test statistic: 109.4194 p-value: 2.276873e-08 Robust Test(5%) : 118.7134 p-value: 9.894441e-10 > ##### Exponentially Weighted Moving Average Approach > da=read.table("m-dec125910-6111.txt",header=T) > head(da) date dec1 dec2 dec5 dec9 dec10 1 19610131 0.058011 0.067392 0.081767 0.096754 0.087207 2 19610228 0.029241 0.042784 0.055524 0.056564 0.060245 3 19610330 0.025896 0.025474 0.041304 0.060563 0.071875 4 19610428 0.005667 0.001365 0.000780 0.011911 0.023328 5 19610531 0.019208 0.036852 0.049590 0.046248 0.050362 6 19610630 -0.024670 -0.025225 -0.040046 -0.050651 -0.051434 > rtn=log(da[,2:4]+1) > m1=VAR(rtn,1) ## Fit VAR(1) model to remove serial correlations Constant term: Estimates: 0.006376978 0.007034631 0.007342962 Std.Error: 0.001759562 0.001950008 0.002237004 AR coefficient matrix AR( 1 )-matrix [,1] [,2] [,3] [1,] -0.194 0.224 0.00836 [2,] -0.232 0.366 -0.04186 [3,] -0.313 0.452 0.00238 standard error [,1] [,2] [,3] [1,] 0.108 0.160 0.101 [2,] 0.120 0.177 0.111 [3,] 0.138 0.204 0.128 Residuals cov-mtx: [,1] [,2] [,3] [1,] 0.001814678 0.001859113 0.001962277 [2,] 0.001859113 0.002228760 0.002420858 [3,] 0.001962277 0.002420858 0.002933081 det(SSE) = 1.712927e-10 AIC = -22.45809 BIC = -22.39289 HQ = -22.43273 > at=m1$residuals ## ARCH test > MarchTest(at) Q(m) of squared series(LM test): Test statistic: 244.7878 p-value: 0 Rank-based Test: Test statistic: 215.215 p-value: 0 Q_k(m) of squared series: Test statistic: 176.9811 p-value: 1.252294e-07 Robust Test(5%) : 155.2633 p-value: 2.347499e-05 > m2=EWMAvol(at,lambda=-0.1) ### Estimation of decaying rate Coefficient(s): Estimate Std. Error t value Pr(>|t|) lambda 0.96427 0.00549 175.6 <2e-16 *** --- > Sigma.t=m2$Sigma.t ### Volatility matrices > m3=MCHdiag(at,Sigma.t) Test results: Q(m) of et: Test and p-value: 59.5436 4.421175e-09 Rank-based test: Test and p-value: 125.0929 0 Qk(m) of epsilon_t: Test and p-value: 189.5403 4.518401e-09 Robust Qk(m): Test and p-value: 228.234 5.57332e-14 > > par(mfcol=c(3,2)) ### Create Figure 7.3 of the text > tdx=c(2:609)/12+1961 > plot(tdx,Sigma.t[,1],xlab='Year',ylab="Variance",type='l') > title(main="(a) Dec 1") > plot(tdx,Sigma.t[,5],xlab='Year',ylab="Variance",type='l') > title(main="(b) Dec 2") > plot(tdx,Sigma.t[,9],xlab='Year',ylab="Variance",type='l') > title(main="(c) Dec 5") > plot(tdx,Sigma.t[,2],xlab='Year',ylab="Covariance",type='l') > title(main="(d) Dec 1 vs Dec 2") > plot(tdx,Sigma.t[,3],xlab='Year',ylab="Covariance",type='l') > title(main="(e) Dec 1 vs Dec 5") > plot(tdx,Sigma.t[,6],xlab='Year',ylab="Covariance",type='l') > title(main="(f) Dec 2 vs Dec 5") > ########## BEKK(1,1) models > da=read.table("m-ibmsp-6111.txt",header=T) > rtn=log(da[,2:3]+1) > m1a=BEKK11(rtn) Initial estimates: 0.00772774 0.005023909 0.06977651 0.02639684 0.03502962 0.1 0.02 0.02 0.1 0.8 0.1 0.1 0.8 Lower limits: -0.0772774 -0.05023909 0.0139553 0.005279367 0.007005923 1e-06 -0.5 -0.5 1e-06 1e-06 -0.5 -0.5 1e-06 Upper limits: 0.0772774 0.05023909 0.07675416 0.02903652 0.03853258 0.999999 0.5 0.5 0.999999 0.999999 0.5 0.5 0.999999 Coefficient(s): Estimate Std. Error t value Pr(>|t|) mu1.ibm 0.00775929 0.00253971 3.05518 0.00224922 ** mu2.sp 0.00565084 0.00154553 3.65624 0.00025594 *** A011 0.01395530 0.00408488 3.41633 0.00063472 *** A021 0.00838972 0.00268086 3.12949 0.00175112 ** A022 0.00700592 0.00193247 3.62537 0.00028855 *** A11 0.15648877 0.06002824 2.60692 0.00913610 ** A21 -0.05926387 0.03253895 -1.82132 0.06855810 . A12 0.23398204 0.08575142 2.72861 0.00636022 ** A22 0.40977179 0.05400961 7.58702 3.2641e-14 *** B11 0.97639151 0.02283328 42.76178 < 2.22e-16 *** B21 0.01449633 0.01196030 1.21204 0.22549813 B12 -0.09287696 0.03227225 -2.87792 0.00400306 ** B22 0.89077633 0.02476925 35.96300 < 2.22e-16 *** --- > names(m1a) [1] "estimates" "HessianMtx" "Sigma.t" > Sigma.t=m1a$Sigma.t > at=cbind(rtn[,1]-0.00776,rtn[,2]-0.00565) > MCHdiag(at,Sigma.t) Test results: Q(m) of et: Test and p-value: 5.280566 0.8716653 Rank-based test: Test and p-value: 16.02931 0.0987965 Qk(m) of epsilon_t: Test and p-value: 29.46281 0.889654 Robust Qk(m): Test and p-value: 52.13744 0.09462895 > ### Figure 7.4 can be drawn using Sigma.t ### ##### Cholesky Decomposition > da=read.table("m-ibmspko-6111.txt",header=T) > rtn=log(da[,2:4]+1) > require(fGarch) > m3=MCholV(rtn) Sample means: 0.007728 0.005024 0.010595 Estimation of the first component Estimate (alpha0, alpha1, beta1): 0.000356 0.117515 0.810288 s.e. : 0.000157 0.037004 0.057991 t-value : 2.262897 3.175772 13.97261 Component 2 Estimation Results (residual series): Estimate (alpha0, alpha1, beta1): 6.4e-05 0.099156 0.858354 s.e. : 3.1e-05 0.027785 0.037238 t-value : 2.034528 3.568616 23.05076 Component 3 Estimation Results (residual series): Estimate (alpha0, alpha1, beta1): 0.000173 0.117506 0.818722 s.e. : 6.2e-05 0.028651 0.038664 t-value : 2.808075 4.101297 21.17521 > names(m3) [1] "betat" "bt" "Vol" "Sigma.t" > at=scale(rtn[37:612,],center=T,scale=F) > Sigma.t=m3$Sigma.t > MCHdiag(at,Sigma.t) Test results: Q(m) of et: Test and p-value: 15.94978 0.1010791 Rank-based test: Test and p-value: 21.99727 0.01511849 Qk(m) of epsilon_t: Test and p-value: 123.7687 0.01057302 Robust Qk(m): Test and p-value: 95.49626 0.3259654 > MCHdiag(at,Sigma.t,5) Test results: Q(m) of et: Test and p-value: 5.717035 0.3347333 Rank-based test: Test and p-value: 5.579834 0.3492711 Qk(m) of epsilon_t: Test and p-value: 59.837 0.06842123 Robust Qk(m): Test and p-value: 58.97874 0.07891472 > #### DCC Models #### Continue to use IBM, SP, and KO log returns #### > m1=dccPre(rtn,include.mean=T,p=0) Sample mean of the returns: 0.00772774 0.005023909 0.01059521 Component: 1 Estimates: 0.000419 0.126739 0.788307 se.coef : 0.000162 0.035405 0.055645 t-value : 2.593448 3.57973 14.16662 Component: 2 Estimates: 9e-05 0.127725 0.836053 se.coef : 4.1e-05 0.03084 0.031723 t-value : 2.20126 4.141592 26.35486 Component: 3 Estimates: 0.000256 0.098705 0.830358 se.coef : 8.5e-05 0.022361 0.033441 t-value : 3.015321 4.414112 24.83088 > names(m1) [1] "marVol" "sresi" "est" "se.coef" > rtn1=m1$sresi > Vol=m1$marVol > m2=dccFit(rtn1) ### Use Tse and Tsui model Estimates: 0.8088046 0.04027251 7.959158 st.errors: 0.1491729 0.02259832 1.135936 t-values: 5.421929 1.782101 7.006697 > names(m2) [1] "estimates" "Hessian" "rho.t" > S2.t=m2$rho.t > m3=dccFit(rtn1,type="Engle") ## Use Engle model Estimates: 0.9126634 0.04530917 8.623668 st.errors: 0.0294762 0.01273911 1.332371 t-values: 30.96272 3.556697 6.472423 > S3.t=m3$rho.t > MCHdiag(rtn1,S2.t) ### Model checking Test results: Q(m) of et: Test and p-value: 20.74259 0.02296172 Rank-based test: Test and p-value: 30.20662 0.0007924436 Qk(m) of epsilon_t: Test and p-value: 132.423 0.002425899 Robust Qk(m): Test and p-value: 109.9671 0.07501565 > MCHdiag(rtn1,S3.t) Test results: Q(m) of et: Test and p-value: 20.02958 0.02897411 Rank-based test: Test and p-value: 27.61638 0.002078829 Qk(m) of epsilon_t: Test and p-value: 131.982 0.002625755 Robust Qk(m): Test and p-value: 111.353 0.06307334 > #### Go-GARCH models > require(gogarch) > crtn=scale(rtn,center=T,scale=F) > m1=gogarch(crtn,~garch(1,1),estby="ica") > m1 **************** *** GO-GARCH *** **************** Components estimated by: fast ICA Dimension of data matrix: (612 x 3). Formula for component GARCH models: ~ garch(1, 1) Orthogonal Matrix U: [,1] [,2] [,3] [1,] 0.3723864 -0.4878111 0.7895370 [2,] -0.7629122 0.3235351 0.5597232 [3,] -0.5284821 -0.8107807 -0.2516770 Linear Map Z: [,1] [,2] [,3] [1,] 0.00873623 -0.03332306 0.060614008 [2,] -0.03042534 -0.00736970 0.030670986 [3,] -0.03975322 -0.04713965 -0.001324943 Estimated GARCH coefficients: omega alpha1 beta1 y1 0.06929475 0.12917711 0.7753224 y2 0.06211986 0.08091875 0.8965807 y3 0.16343143 0.11574017 0.7954333 Convergence codes of component GARCH models: y1 y2 y3 1 1 1 > sigma.t=NULL ### obtain the volatility matrix > for (i in 1:612){ + sigma.t=rbind(sigma.t,c(m1@H[[i]])) + } > MCHdiag(crtn,sigma.t) ### Model checking Test results: Q(m) of et: Test and p-value: 20.99959 0.02109645 Rank-based test: Test and p-value: 30.76079 0.0006425339 Qk(m) of epsilon_t: Test and p-value: 198.4193 3.856975e-10 Robust Qk(m): Test and p-value: 157.3386 1.469088e-05 > M=m1@Z > Minv=solve(M) > bt=crtn%*%t(Minv) > cor(bt^2) [,1] [,2] [,3] [1,] 1.0000000 0.1746403 0.3421938 [2,] 0.1746403 1.0000000 0.1063800 [3,] 0.3421938 0.1063800 1.0000000 > ### Copula approach > da=read.table("m-ibmspko-6111.txt",header=T) > rtn=log(da[,-1]+1) > m1=dccPre(rtn,cond.dist="std") Sample mean of the returns: 0.00772774 0.005023909 0.01059521 Component: 1 Estimates: 0.000388 0.115626 0.805129 9.209269 se.coef : 0.000177 0.036827 0.059471 3.054817 t-value : 2.195398 3.139719 13.5382 3.014671 Component: 2 Estimates: 0.00012 0.130898 0.814531 7.274928 se.coef : 5.7e-05 0.037012 0.046044 1.913331 t-value : 2.102768 3.536655 17.69028 3.802232 Component: 3 Estimates: 0.000216 0.104706 0.837217 7.077138 se.coef : 8.9e-05 0.028107 0.037157 1.847528 t-value : 2.437323 3.725341 22.53208 3.830599 > names(m1) [1] "marVol" "sresi" "est" "se.coef" > Vol=m1$marVol; eta=m1$sresi > > m2=mtCopula(eta,0.8,0.04) Lower limits: 5.1 0.2 1e-04 0.7564334 1.031269 0.8276595 Upper limits: 20 0.95 0.04999999 1.040096 1.417994 1.138032 estimates: 15.38215 0.88189 0.03402486 0.9197241 1.225322 1.058445 std.errors: 8.222771 0.05117132 0.01173272 0.04135667 0.05547635 0.05184866 t-values: 1.870677 17.23407 2.899996 22.23883 22.08729 20.41412 Alternative numerical estimates of se: st.errors: 5.477764 0.05103262 0.01171381 0.04136989 0.05529253 0.05079302 t-values: 2.808107 17.28091 2.904679 22.23173 22.16072 20.83839 > names(m2) [1] "estimates" "Hessian" "rho.t" "theta.t" > MCHdiag(eta,m2$rho.t) Test results: Q(m) of et: Test and p-value: 19.30177 0.03659304 Rank-based test: Test and p-value: 27.03262 0.002573576 Qk(m) of epsilon_t: Test and p-value: 125.9746 0.007387423 Robust Qk(m): Test and p-value: 107.4675 0.1011374 > m3=mtCopula(eta,0.8,0.04,include.th0=F) Value of angles: [1] 0.9455418 1.2890858 1.0345744 Lower limits: 5.1 0.2 1e-05 Upper limits: 20 0.95 0.0499999 estimates: 14.87427 0.8778 0.03365157 std.errors: 7.959968 0.05301289 0.01195093 t-values: 1.868635 16.55824 2.815811 Alternative numerical estimates of se: st.errors: 5.49568 0.05298963 0.01191378 t-values: 2.706539 16.56551 2.824592 > ##### Example 7.6 > da=read.table("d-xomspaapl.txt",header=T) > head(da) date xom sp aapl 1 20070904 0.017497 0.010468 0.041017 2 20070905 -0.000115 -0.011501 -0.051332 3 20070906 0.003096 0.004252 -0.012796 4 20070907 -0.019888 -0.016908 -0.023998 5 20070910 -0.010379 -0.001273 0.037490 6 20070911 0.024511 0.013632 -0.008924 > rtn=log(da[,-1]+1)*100 > mm1=dccPre(rtn,cond.dist="std") Sample mean of the returns: 0.01407631 -0.001785643 0.1231608 Component: 1 Estimates: 0.043949 0.103965 0.882669 6.991655 se.coef : 0.01771 0.020739 0.021429 1.339199 t-value : 2.481573 5.013087 41.19094 5.220775 Component: 2 Estimates: 0.020893 0.104712 0.891949 6.224985 se.coef : 0.009046 0.017481 0.015554 1.24595 t-value : 2.309578 5.990117 57.34556 4.996175 Component: 3 Estimates: 0.064247 0.078172 0.910503 6.702715 se.coef : 0.032488 0.021027 0.023326 1.182064 t-value : 1.977537 3.717595 39.03464 5.670351 > rtn1=mm1$sresi > Vol=mm1$marVol > dim(rtn1) [1] 1280 3 > mm2=mtCopula(rtn1,0.8,0.04) Lower limits: 5.1 0.2 1e-04 0.5675009 0.9203865 0.793023 Upper limits: 20 0.95 0.04999999 0.7803137 1.265531 1.090407 estimates: 8.560638 0.9475189 0.02496858 0.6654369 1.144647 1.027786 std.errors: 1.610043 0.01333286 0.005029271 0.03005827 0.04789454 0.04602053 t-values: 5.317024 71.06644 4.964652 22.13823 23.89932 22.3332 Alternative numerical estimates of se: st.errors: 1.144325 0.0133725 0.005036272 0.02961058 0.04736825 0.04584224 t-values: 7.480947 70.85577 4.957751 22.47294 24.16485 22.42006 > MCHdiag(rtn1,mm2$rho.t) Test results: Q(m) of et: Test and p-value: 24.61657 0.006121662 Rank-based test: Test and p-value: 35.29818 0.0001110462 Qk(m) of epsilon_t: Test and p-value: 171.7599 4.693033e-07 Robust Qk(m): Test and p-value: 113.3854 0.04843521 > #### Principal volatility component analysis > da=read.table("w-7fx.txt",header=T) > head(da) Mon day year GBPUSD NOKUSD SEKUSD CHFUSD CADUSD SGDUSD AUDUSD 1 3 22 2000 0.63397 8.4066 8.6470 1.6560 1.4691 1.7156 1.6480 2 3 29 2000 0.62935 8.4148 8.6232 1.6595 1.4556 1.7152 1.6354 3 4 5 2000 0.62890 8.4670 8.6555 1.6435 1.4543 1.7160 1.6563 4 4 12 2000 0.63029 8.5327 8.6773 1.6450 1.4659 1.7161 1.6744 5 4 19 2000 0.63290 8.6515 8.7601 1.6646 1.4793 1.6999 1.6837 6 4 26 2000 0.63519 8.8426 8.8844 1.7064 1.4763 1.7047 1.6940 > dim(da) [1] 606 10 > fx=log(da[,4:10]) > rtn=diffM(fx) > source("comVol.R") > m1=comVol(rtn,p=5) eigen-values: 4.07592 1.906457 1.646785 1.199888 1.081506 0.7761643 0.4490305 proportion: 0.3660211 0.1712015 0.1478827 0.107751 0.09712014 0.06970022 0.04032333 Checking: Results of individual F-test for ARCH effect Numbers of lags used: 10, 20, 30 Component,(F-ratio P-val) (F-ratio P-val) (F-ratio P-Val) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 21.66 0.00e+00 11.421 0.00e+00 7.725 0.00e+00 [2,] 2 5.78 2.54e-08 4.114 8.66e-09 3.169 7.06e-08 [3,] 3 7.24 7.95e-11 3.804 7.17e-08 2.659 6.66e-06 [4,] 4 4.58 2.76e-06 2.859 3.70e-05 2.279 1.63e-04 [5,] 5 1.49 1.38e-01 0.824 6.85e-01 0.776 7.99e-01 [6,] 6 2.38 9.18e-03 1.872 1.23e-02 1.271 1.55e-01 [7,] 7 1.42 1.68e-01 1.121 3.23e-01 1.126 2.96e-01 > names(m1) [1] "residuals" "values" "vectors" "M" > print(round(m1$M,4)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] -0.1965 -0.2139 -0.0020 0.6560 0.1648 0.3660 -0.2319 [2,] -0.3089 0.0292 -0.1021 0.1588 0.1776 -0.7538 -0.1865 [3,] -0.2351 0.3511 -0.3311 -0.3132 -0.3295 0.5094 -0.2157 [4,] 0.6307 -0.1982 0.1428 -0.0160 -0.2357 -0.0725 -0.2187 [5,] -0.0375 0.0600 0.6978 0.2937 -0.3994 -0.0540 0.5686 [6,] 0.6409 0.8503 -0.0782 0.3649 0.7624 0.1732 0.6629 [7,] 0.0279 -0.2534 0.6055 -0.4761 0.1910 0.0124 -0.2332 > rt=m1$residuals%*%m1$M >