ARCH-GARCH Tutorial with rugarch package
Ozancan Ozdemir
ARCH-GARCH MODELS
The aim of this R tutorial to show when you need (G)ARCH models for volatility and how to fit an appropriate model for your series using rugarch package. Also, you are able to learn how to produce partial bootstrap forecast observations from your GARCH model.
Autoregressive models can be developed for univariate time series data that is stationary (AR), has a trend (ARIMA), and has a seasonal component (SARIMA).One aspect of a univariate time series that these autoregressive models do not model is a change in the variance over time.
Classically, a time series with modest changes in variance can sometimes be adjusted using a power transform, such as by taking the Log or using a Box-Cox transform.There are some time series where the variance changes consistently over time. In the context of a time series in the financial domain, this would be called increasing and decreasing volatility.
In time series where the variance is increasing in a systematic way, such as an increasing trend, this property of the series is called Heteroscedasticity. It’s a fancy word from statistics that means changing or unequal variance across the series.
If the change in variance can be correlated over time, then it can be modeled using an autoregressive or autoregressive moving average process, such as ARCH and GARCH.
What is an ARCH Model?
Autoregressive Conditional Heteroscedasticity, or ARCH, is a method that explicitly models the change in variance over time in a time series.
Specifically, an ARCH method models the variance at a time step as a function of the residual errors from a mean process (e.g. a zero mean).
\(h_t = \omega + \sum_{i}^{q}\alpha_i e^2_{t-i}\)
where \(h_t\) is variance at time t, \(e_{t-i}\) is the model residuals at time t-i.
A lag parameter must be specified to define the number of prior residual errors to include in the model. Using the notation of the GARCH model (discussed later), we can refer to this parameter as “q“.
q: The number of lag squared residual errors to include in the ARCH model.
A generally accepted notation for an ARCH model is to specify the ARCH() function with the q parameter ARCH(q); for example, ARCH(1) would be a first order ARCH model.
The approach expects the series is stationary, other than the change in variance, meaning it does not have a trend or seasonal component. An ARCH model is used to predict the variance at future time steps.
In practice, this can be used to model the expected variance on the residuals after another autoregressive model has been used, such as an ARMA or similar.
“The model should only be applied to a prewhitened residual series \({e_t}\) that is uncorrelated and contains no trends or seasonal changes, such as might be obtained after fitting a satisfactory SARIMA model.” — Page 148, Introductory Time Series with R, 2009.
What is a GARCH Model?
Generalized Autoregressive Conditional Heteroscedasticity, or GARCH, is an extension of the ARCH model that incorporates a moving average component together with the autoregressive component. Specifically, the model includes lag variance terms (e.g. the observations if modeling the white noise residual errors of another process), together with lag residual errors from a mean process.
\(h_t = \omega + \sum_{i}^{q}\alpha_i e^2_{t-i} + \sum_{1}^{p}\beta_i h_{t-i}\)
where \(h_t\) is variance at time t, \(e_{t-i}\) is the model residuals at time t-i.
The introduction of a moving average component allows the model to both model the conditional change in variance over time as well as changes in the time-dependent variance. Examples include conditional increases and decreases in variance.
As such, the model introduces a new parameter “p†that describes the number of lag variance terms:
p: The number of lag variances to include in the GARCH model. q: The number of lag residual errors to include in the GARCH model.
A generally accepted notation for a GARCH model is to specify the GARCH() function with the p and q parameters GARCH(p, q); for example GARCH(1, 1) would be a first order GARCH model.
A GARCH model subsumes ARCH models, where a GARCH(0, q) is equivalent to an ARCH(q) model.
For p = 0 the process reduces to the ARCH(q) process, and for p = q = 0 E(t) is simply white noise. In the ARCH(q) process the conditional variance is specified as a linear function of past sample variances only, whereas the GARCH(p, q) process allows lagged conditional variances to enter as well. This corresponds to some sort of adaptive learning mechanism.
As with ARCH, GARCH predicts the future variance and expects that the series is stationary, other than the change in variance, meaning it does not have a trend or seasonal component.
In the literature, there are lots of GARCH type models.
How to Configure ARCH and GARCH Models
The configuration for an ARCH model is best understood in the context of ACF and PACF plots of the variance of the time series.
This can be achieved by subtracting the mean from each observation in the series and squaring the result, or just squaring the observation if you’re already working with white noise residuals from another model.
The ACF and PACF plots can then be interpreted to estimate values for p and q, in a similar way as is done for the ARMA model.
Example1
Please load the following packages.
library(rugarch)## Loading required package: parallel## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigmalibrary(rmgarch)Remember BDM data set used in the previous recitation. In that lab, we fit \(SARIMA (1,1,1) (0,1,1)_{12}\) model for the series by considering residuals of the model.
bdm=ts(read.table("BDM.txt",header=F),frequency = 12,start=1985)
bdm_train=ts(bdm[1:367],frequency = 12,start=1985)
bdm_test=ts(bdm[368:379],frequency = 12,start=c(2015,8))library(forecast)
fit=Arima(bdm_train,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))
fit## Series: bdm_train
## ARIMA(1,1,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ma1 sma1
## 0.7321 -0.4007 -0.9999
## s.e. 0.1100 0.1584 0.0684
##
## sigma^2 estimated as 0.03884: log likelihood=53.57
## AIC=-99.15 AICc=-99.03 BIC=-83.67r=resid(fit)
#you can also use residuals function to extract residuals of the object.Now, let’s again consider checking heteroscedasticity assumption using visual tools.
ACF & PACF of Squared Residuals
rr=r^2
par(mfrow=c(1,2))
acf(as.vector(rr),main="ACF of Squared Residuals");
pacf(as.vector(rr),main="PACF of Squared Residuals") # homoscedasticity check
Both plots shows that some spikes are out of the white noise bands, i.e there is a correlation between them, that is an indication of heteroscedasticity problem and usage of (G)ARCH model.
However, we need to also check this sitution using a formal test.
Engle’s ARCH Test.
This test is a Lagrange Multiplier test and uses the following hypothesis.
Ho: Residuals exhibits no ARCH effects.
H1: ARCH(lag) effects are present.
library(MTS)
archTest(r)## Q(m) of squared series(LM test):
## Test statistic: 111.6577 p-value: 0
## Rank-based Test:
## Test statistic: 360.9755 p-value: 0Since p values is less than \(\alpha\), we reject Ho. Therefore, we can conclude the presence of ARCH effects.
How to fit model?
There is no certain way to decide best (G)ARCH model for the conditional variance. You will find it by trial and error.
rugarch package
The rugarch package aims to provide for a comprehensive set of methods for modelling univariate GARCH processes, including fitting, filtering, forecasting, simulation as well as diagnostic tools including plots and various tests. Additional methods such as rolling estimation, bootstrap forecasting and simulated parameter density to evaluate model uncertainty provide a rich environment for the modelling of these processes.
In this function, you need to start by specifiying your model. This package seeks and finds the most appropriate GARCH model for your series. To do so, use ugarchspec() function.
Let’s start with default model.
# start with default GARCH spec.
spec = ugarchspec() #the empty function specifies the default model.
print(spec)##
## *---------------------------------*
## * GARCH Model Spec *
## *---------------------------------*
##
## Conditional Variance Dynamics
## ------------------------------------
## GARCH Model : sGARCH(1,1)
## Variance Targeting : FALSE
##
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model : ARFIMA(1,0,1)
## Include Mean : TRUE
## GARCH-in-Mean : FALSE
##
## Conditional Distribution
## ------------------------------------
## Distribution : norm
## Includes Skew : FALSE
## Includes Shape : FALSE
## Includes Lambda : FALSEThis is basic ARIMA(1,1),GARCH(1,1) model.
GARCH(1,1) Model
\(h_t = \omega + \alpha e^2_{t-1} + \beta h_{t-1}\)
library(rugarch)
def.fit = ugarchfit(spec = spec, data = bdm_train)
print(def.fit)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 2.309328 0.130605 17.6818 0.0000
## ar1 0.997292 0.001979 503.8835 0.0000
## ma1 0.450660 0.059689 7.5501 0.0000
## omega 0.000046 0.000041 1.1330 0.2572
## alpha1 0.187315 0.017981 10.4172 0.0000
## beta1 0.811685 0.017057 47.5860 0.0000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 2.309328 0.068412 33.75618 0.000000
## ar1 0.997292 0.004090 243.86034 0.000000
## ma1 0.450660 0.151368 2.97725 0.002908
## omega 0.000046 0.000168 0.27633 0.782291
## alpha1 0.187315 0.040869 4.58328 0.000005
## beta1 0.811685 0.044252 18.34221 0.000000
##
## LogLikelihood : 201.0591
##
## Information Criteria
## ------------------------------------
##
## Akaike -1.06299
## Bayes -0.99914
## Shibata -1.06352
## Hannan-Quinn -1.03762
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 8.862 0.002911
## Lag[2*(p+q)+(p+q)-1][5] 25.518 0.000000
## Lag[4*(p+q)+(p+q)-1][9] 41.119 0.000000
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.003081 0.955734
## Lag[2*(p+q)+(p+q)-1][5] 13.477551 0.001134
## Lag[4*(p+q)+(p+q)-1][9] 15.933786 0.002119
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.3269 0.500 2.000 0.5675
## ARCH Lag[5] 0.4608 1.440 1.667 0.8952
## ARCH Lag[7] 1.5542 2.315 1.543 0.8104
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 6.8626
## Individual Statistics:
## mu 0.2509
## ar1 0.2280
## ma1 0.1576
## omega 0.4691
## alpha1 1.4641
## beta1 0.5673
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.7681 0.07789 *
## Negative Sign Bias 0.6894 0.49101
## Positive Sign Bias 0.1427 0.88659
## Joint Effect 5.8099 0.12124
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 74.91 1.380e-08
## 2 30 87.74 8.073e-08
## 3 40 122.70 1.379e-10
## 4 50 131.77 1.635e-09
##
##
## Elapsed time : 0.3184972In this output, there are several test results related to residuals of the process. Let’s see what these tests do.
First of all, the estimated parameter of ARIMA(1,1) and GARCH(1,1) models were exhibited. It is seen except omega parameter, which is the expected variance value at time t for zero residual and lagged variance values, all parameters are significant.
Then, the some information criterias are given to be used for selecting the best GARCH models for the series. The model having smallest information criteria is the better than the others.
As you would know, Ljung Box Tests are used to test serial autocorrelation among the residuals. (Null: No autocorrelation)
The results show that residuals have autocorrelation, but squared residuals not. (Look p values.)
ARCH LM test is used to check presence of ARCH effect. (Null: Adequately fitted ARCH process)
The results show that the GARCH process is adequately fitted. (Look p values.)
Sign Bias Test is used to test leverage effect in the standardized residuals. (Null: no significiant negative and positive reaction shocks (if exist apARCH type models))
The results show that there is a leverage effect. Fit apARCH type model.
The Nyblom stability test provides a means of testing for structural change within a time series. A structural change implies that the relationship between variables changes overtime e.g. for the regression \(y= \beta x\) beta changes over time. (Null: the parameter values are constant i.e. zero variance, the alternative hypothesis is that their variance > 0.) (Reject Ho, if Test Stat > CV.)
Therefore, we can say that omega, alpha and beta have stability problem. We should also consider TGARCH models.
Adjusted Pearson Goodness-of-Fit Test calculates the chi-squared goodness of fit test, which compares the empirical distribution of the standardized residuals with the theoretical ones from the chosen density. In this case, the chosen density is student t, not standart normal. However, this is not a problem because t-distribution approaches the normal distribution for n>30.
It is seen that we have a normality problem.
Now, by looking at ACF and PACF, we will specify the order of GARCH model manually. The orders are (2,1).
spec=ugarchspec(variance.model = list(model="sGARCH",garchOrder = c(2, 1))) def.fit1= ugarchfit(spec = spec, data = bdm_train)
print(def.fit1)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(2,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 2.296535 0.155024 14.81406 0.000000
## ar1 0.997132 0.002080 479.49953 0.000000
## ma1 0.479947 0.055641 8.62570 0.000000
## omega 0.000075 0.000160 0.46545 0.641610
## alpha1 0.139514 0.038730 3.60217 0.000316
## alpha2 0.066320 0.051571 1.28600 0.198443
## beta1 0.793166 0.039252 20.20686 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 2.296535 0.127004 18.082382 0.000000
## ar1 0.997132 0.006279 158.798243 0.000000
## ma1 0.479947 0.110549 4.341472 0.000014
## omega 0.000075 0.001767 0.042250 0.966300
## alpha1 0.139514 0.585515 0.238276 0.811667
## alpha2 0.066320 0.724459 0.091544 0.927061
## beta1 0.793166 0.487928 1.625581 0.104039
##
## LogLikelihood : 201.4906
##
## Information Criteria
## ------------------------------------
##
## Akaike -1.05989
## Bayes -0.98541
## Shibata -1.06060
## Hannan-Quinn -1.03030
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 7.011 0.0081
## Lag[2*(p+q)+(p+q)-1][5] 23.977 0.0000
## Lag[4*(p+q)+(p+q)-1][9] 39.323 0.0000
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1325 0.715851
## Lag[2*(p+q)+(p+q)-1][8] 12.8869 0.007346
## Lag[4*(p+q)+(p+q)-1][14] 14.7306 0.028283
## d.o.f=3
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[4] 0.08343 0.500 2.000 0.7727
## ARCH Lag[6] 0.23582 1.461 1.711 0.9606
## ARCH Lag[8] 1.91216 2.368 1.583 0.7593
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 7.0422
## Individual Statistics:
## mu 0.2343
## ar1 0.2167
## ma1 0.1623
## omega 0.6559
## alpha1 1.7003
## alpha2 1.0973
## beta1 0.5855
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.80373 0.07211 *
## Negative Sign Bias 0.64647 0.51838
## Positive Sign Bias 0.09539 0.92406
## Joint Effect 6.36512 0.09514 *
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 71.42 5.334e-08
## 2 30 86.76 1.138e-07
## 3 40 112.02 5.549e-09
## 4 50 123.33 2.439e-08
##
##
## Elapsed time : 0.231153Now, we will specify our GARCH model and look for the asymmetric power ARCH model (apARCH) which is the one of the most general type ARCH model. Depending on the values of the parameters, apARCH model can approach another type of (G)ARCH model.
knitr::include_graphics("1.png")
spec=ugarchspec(variance.model = list(model="apARCH"))
#here you can specify also order of the process, and mean model too. spec##
## *---------------------------------*
## * GARCH Model Spec *
## *---------------------------------*
##
## Conditional Variance Dynamics
## ------------------------------------
## GARCH Model : apARCH(1,1)
## Variance Targeting : FALSE
##
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model : ARFIMA(1,0,1)
## Include Mean : TRUE
## GARCH-in-Mean : FALSE
##
## Conditional Distribution
## ------------------------------------
## Distribution : norm
## Includes Skew : FALSE
## Includes Shape : FALSE
## Includes Lambda : FALSEThis is basic ARIMA(1,1),apARCH(1,1) model.
def.fit2= ugarchfit(spec = spec, data = bdm_train)
print(def.fit2)##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : apARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 2.234706 0.031318 71.35483 0.000000
## ar1 0.993135 0.002029 489.52269 0.000000
## ma1 0.441346 0.017122 25.77591 0.000000
## omega 0.021224 0.024897 0.85249 0.393941
## alpha1 0.258540 0.058168 4.44475 0.000009
## beta1 0.771770 0.035488 21.74718 0.000000
## gamma1 -0.613370 0.124872 -4.91198 0.000001
## delta 0.422611 0.289528 1.45965 0.144385
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 2.234706 0.077105 28.98257 0.000000
## ar1 0.993135 0.016405 60.54003 0.000000
## ma1 0.441346 0.111118 3.97186 0.000071
## omega 0.021224 0.194420 0.10917 0.913071
## alpha1 0.258540 0.384820 0.67185 0.501681
## beta1 0.771770 0.173870 4.43878 0.000009
## gamma1 -0.613370 0.828531 -0.74031 0.459112
## delta 0.422611 2.270647 0.18612 0.852351
##
## LogLikelihood : 251.0483
##
## Information Criteria
## ------------------------------------
##
## Akaike -1.3245
## Bayes -1.2394
## Shibata -1.3254
## Hannan-Quinn -1.2907
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 13.43 0.000247
## Lag[2*(p+q)+(p+q)-1][5] 44.18 0.000000
## Lag[4*(p+q)+(p+q)-1][9] 71.57 0.000000
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.6507 0.419880
## Lag[2*(p+q)+(p+q)-1][5] 13.5265 0.001100
## Lag[4*(p+q)+(p+q)-1][9] 16.7278 0.001349
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1834 0.500 2.000 0.6685
## ARCH Lag[5] 0.9849 1.440 1.667 0.7374
## ARCH Lag[7] 2.8863 2.315 1.543 0.5353
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 4.9366
## Individual Statistics:
## mu 2.44123
## ar1 2.68656
## ma1 2.82969
## omega 0.08282
## alpha1 0.08932
## beta1 0.04902
## gamma1 0.15161
## delta 0.36685
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.89 2.11 2.59
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.09146 0.9272
## Negative Sign Bias 1.01795 0.3094
## Positive Sign Bias 0.00168 0.9987
## Joint Effect 1.29982 0.7292
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 63.25 1.174e-06
## 2 30 71.56 1.852e-05
## 3 40 98.29 5.059e-07
## 4 50 95.26 8.404e-05
##
##
## Elapsed time : 0.984658Now, when we look at the information criteria (for e.g AIC or BIC), the apARCH(1,1) model is a better fit than GARCH(1,1) and GARCH(2,1).
After fitting the model, you can obtain useful graphs by using plot function.
plot(def.fit2)Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits 3: Conditional SD (vs |returns|) 4: ACF of Observations 5: ACF of Squared Observations 6: ACF of Absolute Observations 7: Cross Correlation 8: Empirical Density of Standardized Residuals 9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals 11: ACF of Squared Standardized Residuals 12: News-Impact Curve
plot(def.fit2,which=3)
Grey line is the plot of the series. Blue line represents the volatility of the series. By lookng at this plot, you can see the periods in which volatility was high.
Forecasting
There are two forecasting options in this package. Rolling and Bootstrap method.
The bootstrap methods is based on resampling standardized residulas from the emprical distribution of the fitted model to generate future realizations of the series and sigma.
In other words, with bootstrap forecasting, you can forecast both series and conditional variances.
bootp=ugarchboot(def.fit2,method=c("Partial","Full")[1],n.ahead = 12,n.bootpred=1000,n.bootfit=1000)method : Either the full or partial bootstrap (see note).
n.ahead : The forecast horizon.
n.bootfit : The number of simulation based re-fits used to generate the parameter distribution (i.e the parameter uncertainty). Not relevant for the “Partial†method.
n.bootpred: The number of bootstrap replications per parameter distribution per n.ahead forecasts used to generate the predictive density.
bootp##
## *-----------------------------------*
## * GARCH Bootstrap Forecast *
## *-----------------------------------*
## Model : apARCH
## n.ahead : 12
## Bootstrap method: partial
## Date (T[0]): Tem 2015
##
## Series (summary):
## min q.25 mean q.75 max forecast[analytic]
## t+1 6.0733 7.1841 7.3634 7.4921 9.083 7.3625
## t+2 5.2659 7.0081 7.3098 7.5549 11.476 7.3273
## t+3 4.3988 6.8304 7.2509 7.5837 14.862 7.2923
## t+4 3.7313 6.6754 7.1936 7.6100 15.853 7.2576
## t+5 3.1757 6.5462 7.1370 7.6390 13.958 7.2231
## t+6 2.9931 6.4279 7.0764 7.5871 12.163 7.1889
## t+7 2.8320 6.3580 7.0215 7.5500 12.597 7.1549
## t+8 -1.6684 6.2457 6.9749 7.5476 14.280 7.1211
## t+9 -2.9013 6.1955 6.9282 7.4962 16.020 7.0876
## t+10 -1.8986 6.0838 6.8724 7.4404 17.794 7.0542
## .....................
##
## Sigma (summary):
## min q0.25 mean q0.75 max forecast[analytic]
## t+1 0.356818 0.35682 0.35682 0.35682 0.35682 0.35682
## t+2 0.213310 0.28582 0.35059 0.38753 0.81784 0.36470
## t+3 0.155007 0.25217 0.33739 0.38142 1.39115 0.37250
## t+4 0.124663 0.22876 0.33130 0.39462 2.03015 0.38021
## t+5 0.107673 0.21515 0.32721 0.38836 1.50356 0.38782
## t+6 0.072551 0.19736 0.32517 0.39765 1.69172 0.39534
## t+7 0.048986 0.18305 0.31972 0.40074 1.71733 0.40276
## t+8 0.053141 0.17600 0.31627 0.38193 2.28633 0.41008
## t+9 0.045272 0.16646 0.30969 0.38031 2.63000 0.41730
## t+10 0.048916 0.15514 0.30421 0.37559 2.96308 0.42443
## .....................The output above shows the 12 step ahead forecasts for both series and variances.
plot(bootp)Make a plot selection (or 0 to exit):
1: Parameter Density Plots 2: Series Standard Error Plots 3: Sigma Standard Error Plots
plot(bootp,which=2)
Now, let’s compare the performance of apARCH model with SARIMA model fitted at the beginning of the class.
Since the output of bootstrap forecasting methods is S4 class, we need to use “@” sign to take subset.
s_f=bootp@forc@forecast$seriesFor #this is for series forecasts
s_f## Tem 2015
## T+1 7.362513
## T+2 7.327309
## T+3 7.292347
## T+4 7.257625
## T+5 7.223141
## T+6 7.188894
## T+7 7.154881
## T+8 7.121103
## T+9 7.087556
## T+10 7.054240
## T+11 7.021152
## T+12 6.988292s_f1=as.vector(s_f) #for comparison, we make the forecasts as vector.v_f=bootp@forc@forecast$sigmaFor#this is for variance forecasts
v_f## Tem 2015
## T+1 0.3568183
## T+2 0.3647044
## T+3 0.3725011
## T+4 0.3802064
## T+5 0.3878185
## T+6 0.3953357
## T+7 0.4027566
## T+8 0.4100799
## T+9 0.4173045
## T+10 0.4244296
## T+11 0.4314543
## T+12 0.4383781After this, obtain forecast values from the SARIMA model.
f=forecast(fit,h=12)
f## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 2015 7.474373 7.217652 7.731094 7.081752 7.866994
## Sep 2015 7.548222 7.120752 7.975691 6.894464 8.201980
## Oct 2015 7.651434 7.063212 8.239656 6.751827 8.551042
## Nov 2015 7.630804 6.890402 8.371206 6.498456 8.763152
## Dec 2015 7.663498 6.779475 8.547521 6.311501 9.015495
## Jan 2016 7.693069 6.673776 8.712362 6.134195 9.251943
## Feb 2016 7.712925 6.566329 8.859521 5.959358 9.466492
## Mar 2016 7.714172 6.447556 8.980788 5.777050 9.651294
## Apr 2016 7.725560 6.345578 9.105542 5.615060 9.836060
## May 2016 7.719439 6.232123 9.206755 5.444785 9.994092
## Jun 2016 7.672475 6.083268 9.261683 5.241993 10.102958
## Jul 2016 7.725385 6.039189 9.411582 5.146570 10.304200accuracy(f$mean,bdm_test)## ME RMSE MAE MPE MAPE ACF1
## Test set -0.4834464 0.8126902 0.6758118 -7.545091 9.884099 0.5996571
## Theil's U
## Test set 2.029541accuracy(s_f1,bdm_test)## ME RMSE MAE MPE MAPE ACF1
## Test set 0.004245629 0.5087802 0.4164579 -0.5292006 5.815141 0.496337
## Theil's U
## Test set 1.180278As you see, the forecast values obtained from apARCH model outperforms our SARIMA.
REFERENCES
Brownlee, J. (2019, August 21). How to Model Volatility with ARCH and GARCH for Time Series Forecasting in Python. Retrieved from https://machinelearningmastery.com/develop-arch-and-garch-models-for-time-series-forecasting-in-python/.
Ozkan, I. (n.d.) Retrieved from http://yunus.hacettepe.edu.tr/~iozkan/eco665/eco665-17/archgarch.html