# R language uses multivariate ARMA, GARCH, EWMA, ETS, stochastic volatility SV model to model financial time series data

Time：2022-8-5

This article will illustrate different models for univariate and multivariate financial time series, in particular models for conditional means and conditional covariance matrices, volatility.

## mean model

This section explores conditional mean models.

### iid model

We start with a simple iid model. The iid model assumes that the log-return xt is an N-dimensional Gaussian time series:

The sample estimators for the mean and covariance matrix are the sample mean, respectively

and the sample covariance matrix

We start by generating the data, become familiar with the process and ensure that the estimation process gives correct results (i.e. sanity checks). Then use real market data and fit different models.

Let's generate synthetic iid data and estimate mean and covariance matrices:

``# Generate comprehensive income data X &lt;- rmvnorm(n = T, mean = mu, sigma = Sigma)# sample estimates (sample mean and sample covariance matrix) mu\_sm &lt;- colMeans(X)Sigma\_scm &lt;- cov( X)# Error norm(mu\_sm - mu, &quot;2&quot;)#&gt; \[1\] 2.44norm(Sigma\_scm - Sigma, &quot;F&quot;)#&gt; \[1\] 70.79``

Now, let's do it again for a different number of observations T:

``# First generate all data X &lt;- rmvnorm(n = T\_max, mean = mu, sigma = Sigma) # Now iterate over a subset of samples for (T\_ in T\_sweep) { # Sample estimation mu\_sm &lt;- colMeans(X_) Sigma\_scm &lt;- cov(X\_) # Calculate error error\_mu\_vs\_T &lt;- c(error\_mu\_vs\_T, norm(mu\_sm - mu, &quot;2&quot;)) error\_Sigma\_vs\_T &lt;- c(error\_Sigma\_vs\_T, norm(Sigma\_scm - Sigma, &quot;F&quot;))# plot(T\_sweep, error\_mu\_vs\_T, main = &quot;mu estimation error&quot;,``

``plot(T\_sweep, error\_Sigma\_vs\_T main = &quot;Error in Sigma estimation&quot;, ylab = &quot;Error&quot;``

### Univariate ARMA Model

The ARMA(p,q) model on the logarithmic rate of return xt is

where wt is a white noise sequence with zero mean and variance σ2. The parameters of the model are the coefficients ϕi, θi and the noise variance σ2.

Note that the ARIMA(p,d,q) model is an ARMA(p,q) model with a time difference of order d. So if we replace the log price with xt, the previous log return model is actually an ARIMA(p, 1, q) model because once the log price is differentiated, we get the log return.

#### rugarchGenerate data

we will userugarchThe package generates univariate ARMA data, estimates parameters and makes predictions.

First, we need to define the model:

``# Specify AR(1) model with given coefficients and parameters #&gt; #&gt; *------------------------------- ---*#&gt; * ARFIMA Model Spec *#&gt; *------------------------------------------------*# &gt; Conditional Mean Dynamics#&gt; ------------------------------------#&gt; Mean Model : ARFIMA(1 ,0,0)#&gt; Include Mean : TRUE #&gt; #&gt; Conditional Distribution#&gt; ----------------------------------------- -----#&gt; Distribution : norm #&gt; Includes Skew : FALSE #&gt; Includes Shape : FALSE #&gt; Includes Lambda : FALSE#&gt; Level Fixed Include Estimate LB UB#&gt; mu 0.01 1 1 0 NA NA#&gt; ar1 -0.90 1 1 0 NA NA#&gt; ma 0.00 0 0 0 NA NA#&gt; arfima 0.00 0 0 0 NA NA#&gt; archm 0.00 0 0 0 NA NA#&gt; mxreg 0.00 0 0 0 NA NA#&gt; sigma 0.20 1 1 0 NA NA #&gt; alpha 0.00 0 0 0 NA NA#&gt; beta 0.00 0 0 0 NA NA#&gt; gamma 0.00 0 0 0 NA NA#&gt; eta1 0.00 0 0 0 NA NA#&gt; eta2 0.00 0 0 0 NA NA#&gt; delta 0.00 0 0 0 NA NA#&gt; lambda 0.00 0 0 0 NA NA#&gt; vxreg 0.00 0 0 0 NA NA#&gt; skew 0.00 0 0 0 NA NA#&gt; shape 0.00 0 0 0 NA NA#&gt; ghlambda 0.00 0 0 0 NA NA#&gt; xi 0.00 0 0 0 NA NA NAfixed.pars#&gt; \$mu#&gt; \[1\] 0.01#&gt; #&gt; \$ar1#&gt; \[1\ ] -0.9#&gt; #&gt; \$sigma#&gt; \[1\] 0.2true_params#&gt; mu ar1 sigma #&gt; 0.01 -0.90 0.20``

Then, we can generate the time series:

``# Simulate a path apath(spec, n.sim = T)# Convert to xts and plot plot(synth\_log\_returns, main = &quot;Log returns for ARMA model&quot; plot(synth\_log\_prices, main = &quot; Logarithmic price of ARMA model&quot;``

#### ARMA model

Now, we can estimate the parameters (which we already know):

``# Specify AR(1) model arfimaspec(mean.model = list(armaOrder = c(1,0), include.mean = TRUE))# Estimation model #&gt; mu ar1 sigma #&gt; 0.0083 -0.8887 0.1987#&gt; mu ar1 sigma #&gt; 0.01 -0.90 0.20``

We can also study the effect of the sample size T on the parameter estimation error:

``# loop for (T_ in T\_sweep) { estim\_coeffs\_vs\_T &lt;- rbind(estim\_coeffs\_vs\_T, coef(arma\_fit)) error\_coeffs\_vs\_T &lt;- rbind(error\ _coeffs\_vs\_T, abs(coef(arma\_fit) - true\_params)/true\_params)# plot matplot(T\_sweep, estim\_coeffs\_vs_T, main = &quot;Estimated ARMA coefficients&quot;, xlab = &quot; T&quot;, ylab = &quot;value&quot;,``

``matplot(T\_sweep, 100*error\_coeffs\_vs\_T, main = &quot;Relative error in estimating ARMA coefficients&quot;, xlab = &quot;T&quot;, ylab = &quot;Error (%)&quot;,``

First, the true μ is almost zero, so the relative error may appear unstable. After T = 800 samples, the other coefficients are well estimated.

#### ARMApredict

For a sanity check, we will now compare the two packagesForecastandrugarch results

``# specify AR(1) model with given coefficients and parameters spec(mean.model = list(armaOrder = c(1,0), include.mean = TRUE), fixed.pars = list(mu = 0.005, ar1 = -0.9, sigma = 0.1))# Generate a sequence of length 1000 arfima(arma\_fixed\_spec, n.sim = 1000)@path\$seriesSim # Specify and fit a model using the rugarch package spec(mean.model = list( armaOrder = c(1,0), include.mean = TRUE))# Fit a model using package &quot;forecast&quot; #&gt; ARIMA(1,0,0) with non-zero mean #&gt; #&gt; Coefficients:#&gt; ar1 mean #&gt; -0.8982 0.0036#&gt; se 0.0139 0.0017#&gt; #&gt; sigma^2 estimated as 0.01004: log likelihood=881.6#&gt; AIC=-1757.2 AICc=-1757.17 BIC=-1742.47# Compare model coefficients#&gt; ar1 intercept sigma #&gt; -0.898181148 0.003574781 0.100222964#&gt; mu ar1 sigma #&gt; 0.003605805 -0.898750138 0.100199956``

Indeed, both packages gave the same results.

#### ARMA model selection

In previous experiments, we assumed that we knew the order of the ARMA model, i.e. p = 1 and q = 0. In reality, the order is unknown, so different order combinations have to be tried. The higher the order, the better the fit, but this will inevitably lead to overfitting. Many methods have been developed to penalize the increase in complexity to avoid overfitting, such as AIC, BIC, SIC, HQIC, etc.

``# Try different combinations# View ranking#&gt; AR MA Mean ARFIMA BIC converged#&gt; 1 1 0 1 0 -0.38249098 1#&gt; 2 1 1 1 0 -0.37883157 1#&gt; 3 2 0 1 0 -0.37736340 1#&gt; 4 1 2 1 0 -0.37503980 1#&gt; 5 2 1 1 0 -0.37459177 1#&gt; 6 3 0 1 0 -0.37164609 1#&gt; 7 1 3 1 0 -0.37143480 1#&gt; 8 2 2 1 0 -0.37107841 1#&gt; 9 3 1 1 0 -0.36795491 1#&gt; 10 2 3 1 0 -0.36732669 1#&gt; 11 3 2 1 0 -0.36379209 1#&gt; 12 3 3 1 0 -0.36058264 1#&gt; 13 0 3 1 0 -0.11875575 1#&gt; 2 1 0 0.02957266 1#&gt; 15 0 1 1 0 0.39326050 1#&gt; 16 0 0 1 0 1.17294875 1# choose the best armaOrder#&gt; AR MA #&gt; 1 0``

In this case, since the number of observations T = 1000 is large enough, the order is correctly detected. Conversely, if you try to use T=200, the detected order is p=1, q=3.

#### ARMA forecast

Once the ARMA model parameters ϕi^i and θ^j have been estimated, the model can be used to predict future values. For example, the prediction of xt based on past information is

And the prediction error will be xt-x^t=wt (assuming the parameters have been estimated) with a variance of σ2. packagerugarchMakes predictions on out-of-sample data simple:

``# Estimated model (excluding out-of-sample) coef(arma\_fit)#&gt; mu ar1 sigma #&gt; 0.007212069 -0.898745183 0.200400119# Predicted log-returns out-of-sample forecast\_log\_returns &lt;- xts(arma\[email protected]\$ seriesFor\[1, \], dates\_out\_of\_sample)# restore log price prev\_log\_price &lt;- head(tail(synth\_log\_prices, out\_of\_sample+1), out\_of \_sample)# log-return graph plot(cbind(&quot;fitted&quot; = fitted(arma\_fit),# log-price graph plot(cbind(&quot;forecast&quot; = forecast\_log\_prices, main = &quot;log price forecast&quot; , legend.loc = &quot;topleft&quot;)``

### Multivariate VARMA Model

The VARMA(p,q) model on the logarithmic rate of return xt is

where wt is a white noise sequence with zero mean and covariance matrix Σw. The parameters of this model are the vector/matrix coefficients ϕ0, Φi, Θj and the noise covariance matrix Σw.

### Compare

``# Load S&amp;P 500 data head(SP500\_index\_prices)#&gt; SP500#&gt; 2012-01-03 1277.06#&gt; 2012-01-04 1277.30#&gt; 2012-01-05 1281.06#&gt; 2012-01-06 1277.81# &gt; 2012-01-09 1280.70#&gt; 2012-01-10 1292.08# Prepare training and test data logreturns\_trn &lt;- logreturns\[1:T\_trn\]logreturns\_tst &lt;- logreturns\[-c(1:T \_trn)\]# Plot { plot(logreturns, addEventLines(xts(&quot;training&quot;``

Now, we use the training data (i.e., for t=1,…,Ttrnt=1,…,Ttrn) to fit different models (note that out-of-sample data is excluded by indicating`out.sample = T_tst`). In particular, we will consider iid models, AR models, ARMA models, and some ARCH and GARCH models (variance modeling will be studied in more detail later).

``# Fit the iid model coef(iid\_fit)#&gt; mu sigma #&gt; 0.0005712982 0.0073516993mean(logreturns\_trn)#&gt; \[1\] 0.0005681388sd(logreturns\_trn)#&gt; \[1\] 0.007360208# Fitting AR(1) model coef(ar\_fit)#&gt; mu ar1 sigma #&gt; 0.0005678014 -0.0220185181 0.0073532716# Fit ARMA(2,2) model coef(arma\_fit)#&gt; mu ar1 ar2 ma1 ma2 sigma #&gt; 0.0007223304 0.0268612 0.9095552008 -0.0832923604 -0.9328475211 0.0072573570# Fit ARMA(1,1) + ARCH(1) model coef(arch\_fit)#&gt; mu ar1 ma1 omega alpha1 #&gt; 6.321441e-04 8.720929e-02 -9.398810 e-05 9.986975e-02#Fit ARMA(0,0)+ARCH(10) model coef(long\_arch\_fit)#&gt; mu omega alpha1 alpha2 alpha3 alpha4 alpha5 #&gt; 7.490786e-04 2.452099e-05 6.888561 e-02 7.207551e-02 1.419938e-01 1.909541e-02 3.082806e-02 #&gt; alpha6 alpha7 alpha8 alpha9 alpha10 #&gt; 4.026539e-02 3.050040e-07 9.260183e -02 1.150128e-01 1.068426e-06# Fit ARMA(1,1)+GARCH(1,1) model coef(garch_fit)#&gt; mu ar1 ma1 omega alpha1 beta1 #&gt; 6.660346e-04 9.664597e-01 - 1.000000e+00 7.066506e-06 1.257786e-01 7.470725e-01``

We use different models to predict log returns:

``# Prepare to forecast logarithmic returns for out-of-sample periods # iid model forecast forecast(iid\_fit, n.ahead = 1, n.roll = T\_tst - 1) dates\_out\_of\_sample)# AR(1) model Make prediction forecast(ar\_fit, n.ahead = 1, n.roll = T\_tst - 1) dates\_out\_of\_sample)# ARMA(2,2) model make prediction forecast(arma\_fit, n. ahead = 1, n.roll = T\_tst - 1) dates\_out\_of\_sample)# Use ARMA(1,1) + ARCH(1) model for prediction forecast(arch\_fit, n.ahead = 1, n.roll = T\_tst - 1) dates\_out\_of\_sample)# ARMA(0,0)+ARCH(10) model forecast forecast(long\_arch\_fit, n.ahead = 1, n.roll = T\_tst - 1) dates\_out\_of\_sample)# ARMA(1,1)+GARCH(1,1) model forecast forecast(garch\_fit, n.ahead = 1, n.roll = T\_tst - 1) dates\_out\_of_sample)``

We can calculate the prediction errors (in-sample and out-of-sample) for different models:

``print(error_var)#>                           in-sample out-of-sample#> iid                    5.417266e-05  8.975710e-05#> AR(1)                  5.414645e-05  9.006139e-05#> ARMA(2,2)              5.265204e-05  1.353213e-04#> ARMA(1,1) + ARCH(1)    5.415836e-05  8.983266e-05#> ARCH(10)               5.417266e-05  8.975710e-05#> ARMA(1,1) + GARCH(1,1) 5.339071e-05  9.244012e-05``

We can observe that as the model complexity increases, the in-sample error tends to be smaller (due to higher degrees of freedom in fitting the data), although the difference is negligible. What matters is actually the out-of-sample error: we can see that increasing the model complexity may give poorer results. In terms of error in predicting returns, it seems that the simplest iid model is sufficient.

Finally, let's show some graphs of out-of-sample error:

``plot(error, main = &quot;Out-of-Sample Error in Return Forecasts for Different Models&quot;,``

Note that since we did not refit the model, the error increases over time (especially for ARCH modeling).

### Scrolling window comparison

Let's first compare the concepts of static forecasting versus rolling forecasting with a simple example:

``#ARMA(2,2) model spec &lt;- spec(mean.model = list(armaOrder = c(2,2), include.mean = TRUE))# static fit and prediction ar\_static\_fit &lt;- fit( spec = spec, data = logreturns, out.sample = T\_tst)#Roll fit and prediction modelroll &lt;- aroll(spec = spec, data = logreturns, n.ahead = 1, #Prediction plot plot(cbind(&quot;static forecast&quot; = ar\_static\_fore\_logreturns, main = &quot;Prediction using ARMA(2,2) model&quot;, legend.loc = &quot;topleft&quot;)# forecast error plot plot(error_logreturns, col = c(&quot;black&quot;, &quot;red&quot;), lwd = 2, main = &quot;Prediction error of ARMA(2,2) model&quot;, legend.loc = &quot;topleft&quot;)``

We can clearly observe the effect of the rolling window process on the time series.

Now, we can redo all predictions for all models on a rolling window basis:

``# Rolling forecast based on iid model roll(iid\_spec, data = logreturns, n.ahead = 1, forecast.length = T\_t# Rolling forecast of AR(1) model roll(ar\_spec, data = logreturns, n .ahead = 1, forecast.length = T\_tst, # Rolling forecast for ARMA(2,2) model roll(arma\_spec, data = logreturns, n.ahead = 1, forecast.length = T\_tst, # ARMA (1,1) + ARCH(1) model rolling forecast roll(arch\_spec, data = logreturns, n.ahead = 1, forecast.length = T\_tst, refit.every = 50, refit.win# ARMA( 0,0) + ARCH(10) model rolling forecast roll(long\_arch\_spec, data = logreturns, n.ahead = 1, forecast.length = T\_tst, refit.every = 50, # ARMA(1, 1) + GARCH(1,1) model rolling forecast roll(garch\_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst, refit.every = 50, refit.window``

Let's look at the forecast error in the rolling benchmark case:

``print(rolling\_error\_var)#>                           in-sample out-of-sample#> iid                    5.417266e-05  8.974166e-05#> AR(1)                  5.414645e-05  9.038057e-05#> ARMA(2,2)              5.265204e-05  8.924223e-05#> ARMA(1,1) + ARCH(1)    5.415836e-05  8.991902e-05#> ARCH(10)               5.417266e-05  8.976736e-05#> ARMA(1,1) + GARCH(1,1) 5.339071e-05  8.895682e-05``

and some graphs:

``plot(error_logreturns, main = &quot;Rolling prediction errors for different models&quot;, legend.loc = &quot;topleft&quot;``

We see that now all models fit the time series. Furthermore, we did not find any significant differences between the models.

We can finally compare static and rolling errors:

``barplot(rbind(error\_var\[, &quot;out-of-sample&quot;\], rolling\_error_var\[, &quot;out-of-sample&quot;\]) col = c(&quot;darkblue&quot;, &quot;darkgoldenrod&quot;), legend = c(&quot;Static Forecast&quot;, &quot;Rolling Forecast&quot;),``

We can see that rolling forecasts are necessary in some cases. So, in practice, we need rolling forecast improvements on a regular basis.

## variance model

### ARCH and GARCH models

The ARCH(m) model for the log-return residual wt is

where zt is a white noise sequence with zero mean and constant variance, and the conditional variance σ2t is modeled as

Among them, m is the model order, ω&gt;0, αi≥0 are parameters.

The GARCH(m,s) model extends the ARCH model with a recursive term on σ2t:

The parameters ω\&gt; 0, αi≥0, βj≥0 need to satisfy the stability of ∑mi = 1αi+ ∑sj = 1βj≤1.

#### rugarchGenerate data

First, we need to define the model:

``# Specify a GARCH model with given coefficients and parameters #&gt; #&gt; *---------------------------------* #&gt; * GARCH Model Spec *#&gt; *---------------------------------*#&gt; #&gt; Conditional Variance Dynamics #&gt; ------------------------------------#&gt; GARCH Model : sGARCH(1,1) #&gt; Variance Targeting : FALSE #&gt; #&gt; Conditional Mean Dynamics#&gt; ----------------------------------- -#&gt; Mean Model : ARFIMA(1,0,0)#&gt; Include Mean : TRUE #&gt; GARCH-in-Mean : FALSE #&gt; #&gt; Conditional Distribution#&gt; ------------- -----------------------#&gt; Distribution : norm #&gt; Includes Skew : FALSE #&gt; Includes Shape : FALSE #&gt; Includes Lambda : FALSE#&gt; Level Fixed Include Estimate LB UB#&gt; mu 0.005 1 1 0 NA NA#&gt; ar1 -0.900 1 1 0 NA NA#&gt; ma 0.000 0 0 0 NA NA#&gt; arfima 0.000 0 0 0 NA NA#&gt; archm 0.000 0 0 0 NA NA #&gt; mxreg 0.000 0 0 0 NA NA#&gt; omega 0.001 1 1 0 NA NA #&gt; alpha1 0.300 1 1 0 NA NA#&gt; beta1 0.650 1 1 0 NA NA#&gt; gamma 0.000 0 0 0 NA NA#&gt; eta1 0.000 0 0 0 NA NA#&gt; eta2 0.000 0 0 0 NA NA#&gt; delta 0.000 0 0 0 NA NA#&gt; lambda 0.000 0 0 0 NA NA#&gt; vxreg 0.000 0 0 0 NA NA#&gt; skew 0.000 0 0 0 NA NA#&gt; shape 0.000 0 0 0 NA NA#&gt; ghlambda 0.000 0 0 0 NA NA# &gt; xi 0.000 0 0 0 NA NA#&gt; \$mu#&gt; \[1\] 0.005#&gt; #&gt; \$ar1#&gt; \[1\] -0.9#&gt; #&gt; \$omega#&gt; \[1\] 0.001# &gt; #&gt; \$alpha1#&gt; \[1\] 0.3#&gt; #&gt; \$beta1#&gt; \[1\] 0.65true_params#&gt; mu ar1 omega alpha1 beta1 #&gt; 0.005 -0.900 0.001 0.300 0.650``

Then, we can generate the returns time series:

``# Simulate a path hpath(garch\_spec, n.sim = T)#&gt; num \[1:2000, 1\] 0.167 -0.217 # Plot log returns { plot(synth\_log\_returns, main = &quot;GARCH model &quot;, lwd = 1.5) lines(synth\_volatility``

#### GARCH

Now, we can estimate the parameters:

``# Specify a GARCH model ugarchspec(mean.model = list(armaOrder = c(1,0)# Estimate the model coef(garch_fit)#&gt; mu ar1 omega alpha1 beta1 #&gt; 0.0036510100 -0.8902333595 0.00081811434 0.2810460728 0.6717486402 omega ar beta1 o #&gt; 0.005 -0.900 0.001 0.300 0.650# Coefficient error#&gt; mu ar1 omega alpha1 beta1 #&gt; 0.0013489900 0.0097666405 0.0001188566 0.0189539272 0.0217486402``

We can also study the effect of the sample size T on the parameter estimation error:

``# 循环for (T_ in T\_sweep) {  garch\_fit   error\_coeffs\_vs\_T <- rbind(error\_coeffs\_vs\_T, abs((coef(garch\_fit) - true\_params)/true\_params))  estim\_coeffs\_vs\_T <- rbind(estim\_coeffs\_vs\_T, coef(garch\_fit))# 绘图matplot(T\_sweep, 100*error\_coeffs\_vs\_T,         main = "估计GARCH系数的相对误差", xlab = "T", ylab = "误差 (%)",``

The true ω is almost zero, so the error is very unstable. As for the other coefficients, as in the ARMA case, the estimates of μ are really poor (more than 50% relative error), while the other coefficients seem to be well estimated after T=800 samples.

#### Comparison of GARCH results

As a sanity check, we will now compare the two packagesfGarchandrugarch results

``# Specify ARMA(0,0)-GARCH(1,1) with specific parameter values ​​as data generation process garch\_spec #Generate data with length 1000 path(garch\_fixed\_spec, n.sim = 1000)@path \$# Specify and fit model using 'rugarch' package rugarch\_fit &lt;- ugarchfit(spec = garch\_spec, data = x) # Fit model using 'fGarch' package garchFit(formula = ~ garch(1, 1), data = x, trace = FALSE)# compare model coefficients #&gt; mu omega alpha1 beta1 #&gt; 0.09749904 0.01395109 0.13510445 0.73938595#&gt; mu omega alpha1 beta1 #&gt; 0.09750394 0.0139G2648 0.13527024 compare standard deviation of fit f5_839716 #&gt; \[1\] 0.3513549 0.3254788 0.3037747 0.2869034 0.2735266 0.2708994``

Indeed, both packages gave the same results.

#### userugarchpackage for GARCH prediction

Once the parameters of a GARCH model have been estimated, the model can be used to predict future values. For example, a one-step prediction of conditional variance based on past information is

Given ω^/(1-∑mi=1α^i-∑sj=1β^j). packagerugarchMakes predictions on out-of-sample data simple:

``# Estimate the model, excluding out-of-sample garch\_fit coef(garch\_fit)#&gt; mu ar1 omega alpha1 beta1 #&gt; 0.0034964331 -0.8996287630 0.0006531088 0.3058756796 0.6815452241# Predict the log return for the entire sample garch\[email protected]\$sigmaFor\[forecast\$sigmaFor\ , \]# log-return plot plot(cbind(&quot;fitted&quot; = fitted(garch\_fit), main = &quot;synthetic log-return forecast&quot;, legend.loc = &quot;topleft&quot;)``

``#Volatility log-return plot plot(cbind(&quot;fitted volatility&quot; = sigma(garch_fit), main = &quot;Predicted volatility of synthetic log-return&quot;, legend.loc = &quot;topleft&quot;)``

### different methods

``# Load S&amp;P 500 index data head(SP500\_index\_prices)#&gt; SP500#&gt; 2008-01-02 1447.16#&gt; 2008-01-03 1447.16#&gt; 2008-01-04 1411.63#&gt; 2008-01-07 1416.18 #&gt; 2008-01-08 1390.19#&gt; 2008-01-09 1409.13# Prepare training and test data x\_trn &lt;- x\[1:T\_trn\]x\_tst &lt;- x\[-c(1: T\_trn)\]# plot { plot(x, main = &quot;revenue&quot; addEventLines(xts(&quot;training&quot;, in``

#### constant

``plot(cbind(sqrt(var\_constant), x\_trn) main = &quot;constant&quot;)``

#### moving average

Now, let's use a moving average of squared returns:

``plot(cbind(sqrt(var\_t), x\_trn), main = &quot;Envelope based on simple rolling square mean (time period=20)``

#### EWMA

Exponentially Weighted Moving Average (EWMA):

Note that this can also be modeled as an ETS(A,N,N) state space model:

``plot(cbind(std\_t, x\_trn), main = &quot;Envelope based on squared EWMA&quot;)``

#### Multiplication ETS

We can also try different variants of the ETS model. For example, a multiplicative noise version ETS(M,N,N) with a state-space model:

``plot(cbind(std\_t, x\_trn), col = c(&quot;red&quot;, &quot;black&quot;) main = &quot;Envelope of square-based ETS(M,N,N)&quot;``

#### ARCH

Now, we can use more complex ARCH modeling:

``plot(cbind(std\_t, x\_trn), col = c(&quot;red&quot;, &quot;black&quot;) main = &quot;ARCH(5) based envelope&quot;)``

#### GARCH

We can lift the model to GARCH:

``plot(cbind(std\_t, x\_trn), col = c(&quot;red&quot;, &quot;black&quot;) main = &quot;GARCH(1,1) based envelope&quot;)``

#### SV Stochastic Volatility

Finally, we can use stochastic volatility modeling:

or, equivalently,

``plot(cbind(std\_t, x\_trn), col = c(&quot;red&quot;, &quot;black&quot;), main = &quot;Envelope Analysis Based on Stochastic Volatility&quot;)``

#### Compare

We can now compare the error in the variance estimates for each method over the out-of-sample period:

``#&gt; MA EWMA ETS(M,N,N) ARCH(5) GARCH(1,1) SV #&gt; 2.204965e-05 7.226188e-06 3.284057e-06 7.879039e-05 6.496545e-06 6.705059e-06barplot( error_all, main = &quot;Error in out-of-sample variance estimation&quot;``

### Scrolling window comparison

Rolling window comparison of six methods: MA, EWMA, ETS(MNN), ARCH(5), GARCH(1,1) and SV.

``#rolling window lookback &lt;- 200len\_tst &lt;- 40for (i in seq(lookback, T-len\_tst, by = len\_tst)) { # MA var\_t &lt;- roll\_meanr(x\_trn^2, n = 20, fill = NA) var\_fore &lt;- var(x\_trn/sqrt(var\_t), na.rm = TRUE) * tail(var\_t, 1) error\_ma &lt;- c(error\ _ma, abs(var\_fore - var\_tst)) # EWMA error\_ewma &lt;- c(error\_ewma, abs(var\_fore - var\_tst)) # ETS(M,N,N) error\_ets\ _mnn &lt;- c(error\_ets\_mnn, abs(var\_fore - var\_tst)) # ARCH error\_arch &lt;- c(error\_arch, abs(var\_fore - var\_tst)) # GARCH error\ _garch &lt;- c(error\_garch, abs(var\_fore - var\_tst)) # SV error\_sv &lt;- c(error\_sv, abs(var\_fore - var\_tst))}barplot(error_all, main = &quot;Variance Estimation Error&quot;,``

### Multivariate GARCH Model

For illustration purposes, we will only consider the Constant Conditioned Correlation (CCC) and Dynamic Conditioned Correlation (DCC) models as they are the most popular. The logarithmic rate of return residual wt is modeled as

where zt is an iid white noise sequence with zero mean and constant covariance matrix II. The conditional covariance matrix Σt is modeled as

where Dt = Diag(σ1,t,…,σN,t) is the normalized noise vector C, and the covariance matrix ηt=C-1wt (ie, it contains diagonal elements equal to 1).

Basically, with this model, the diagonal matrix Dt contains a set of univariate GARCH models, and then the matrix C contains some correlations between the sequences. The main disadvantage of this model is that the matrix C is constant. To overcome this problem, DCC is proposed as

where Ct contains diagonal elements equal to 1. To enforce a diagonal element equal to 1, Engle models it as

Qt has arbitrary diagonal elements and follows the model

We will generate data, estimate parameters and forecasts.

• SPDR S＆P 500 ETF
• 20+ Year Treasury Bond ETFs
• IEF: 7-10 Year Treasury Bond ETF
``# Download data prices &lt;- xts()head(prices)#&gt; SPY TLT IEF#&gt; 2013-01-02 127.8779 99.85183 93.65224#&gt; 2013-01-03 127.5890 98.49886 93.17085#&gt; 2013-01-04 126# 9&gt;3.849 2013-01-07 127.7991 98.92480 93.26714#&gt; 2013-01-08 127.4314 99.57622 93.49468#&gt; 2013-01-09 127.7553 99.48438 93.54719# Plot three &quot;logarithmic price series of three ETFs&quot; Number price&quot;, legend.loc = &quot;topleft&quot;)``

First, we define the model:

``# Specify iid univariate time series model ugarch_spec # Specify DCC model spec( multispec(replicate(spec, n = 3))``

Next, we fit the model:

``# 估计模型#> #> *---------------------------------*#> *          DCC GARCH Fit          *#> *---------------------------------*#> #> Distribution         :  mvnorm#> Model                :  DCC(1,1)#> No. Parameters       :  44#> \[VAR GARCH DCC UncQ\] : \[30+9+2+3\]#> No. Series           :  3#> No. Obs.             :  1007#> Log-Likelihood       :  12198.4#> Av.Log-Likelihood    :  12.11 #> #> Optimal Parameters#> -----------------------------------#>               Estimate  Std. Error   t value Pr(>|t|)#> \[SPY\].omega   0.000004    0.000000  11.71585 0.000000#> \[SPY\].alpha1  0.050124    0.005307   9.44472 0.000000#> \[SPY\].beta1   0.870051    0.011160  77.96041 0.000000#> \[TLT\].omega   0.000001    0.000001   0.93156 0.351563#> \[TLT\].alpha1  0.019716    0.010126   1.94707 0.051527#> \[TLT\].beta1   0.963760    0.006434 149.79210 0.000000#> \[IEF\].omega   0.000000    0.000001   0.46913 0.638979#> \[IEF\].alpha1  0.031741    0.023152   1.37097 0.170385#> \[IEF\].beta1   0.937777    0.016498  56.84336 0.000000#> \[Joint\]dcca1  0.033573    0.014918   2.25044 0.024421#> \[Joint\]dccb1  0.859787    0.079589  10.80278 0.000000#> #> Information Criteria#> ---------------------#>                     #> Akaike       -24.140#> Bayes        -23.925#> Shibata      -24.143#> Hannan-Quinn -24.058#> #> #> Elapsed time : 0.8804049``

We can plot the time-varying correlation:

``# Extract time-varying covariance and correlation matrix dim(dcc\_cor)#&gt; \[1\] 3 3 1007#plot(corr\_t main = &quot;time-varying correlation&quot;, legend.loc = &quot;left&quot;)``

We see a very high and fairly stable correlation between the two income ETFs. The correlation with SPY is small and fluctuates in an interval less than 0.

Most Popular Insights

2.HAR-RV Model Based on Mixed Data Sampling (MIDAS) Regression in R Language to Predict GDP GrowthThe regression HAR-RV model predicts GDP growth&quot;)

## Map

1. Map Map is an unordered key-value based data structure. Map in GO language is a reference type and must be initialized with Make before it can be used. 1.1.1 map definition The definition syntax of map in GO language is as follows map[KeyType]ValueType in, KeyType: Indicates the type of key. ValueType: Indicates the type […]