Thursday, July 23, 2015

GARCH model estimation, Backtesting the risk model and Forecasting

In this blog I will show you the GARCH model estimation, Backtesting the risk model and Forecasting. Here is my Video:


STEP-1: Installing required packages and loading them as library in R

We will need following packages so if you have not installed them use the following command.
install.packages(c("tseries", "zoo", "forecast", "FinTS", "rugarch"))

Now, load these packages in your R.
library("tseries")
library("zoo")
library("forecast")
library("FinTS")
library("rugarch")

STEP-2: Downloading the Data

Let’s download the daily adjusted close data of Dow Jones Industrial Average (^DJI) from first of January 2010 till 23rd of July 2015 (today is 24th July) using Yahoo then return that to zoo class object named DJI.data.

DJI.data = get.hist.quote(instrument="^DJI", start="2010-01-01",  end="2015-07-23", quote="AdjClose", provider="yahoo", compression="d", retclass="zoo")

You can use following command to plot that DJI.data object.
plot(DJI.data, main = "DJIE Closing Prices on NASDAQ", ylab = "Price (USD)", xlab = "Date")

Fig-1

STEP-3: Getting the Returns of the Data

Let’s take the log return of DJI.data and name the object ans DJI.ret and plot the data. You can see in the Fig-2 that high volatile days are followed by high volatile days and low volatile days by low volatile days.

DJI.ret <- diff(log(DJI.data)) * 100
plot(DJI.ret, main = "Monthly Compound Returns", xlab = "Date", ylab = "Return in percent")

Fig-2

STEP-4: Finding Best Mean Model Using ARIMA

Now, let’s try to fit a best ARIMA (p,d,q) model based upon the Bayesian Information Criterions. Use the following command. Note that returns are always reverting to some mean returns hence you should expect the returns to be I(0). Hence, the ARIMA is actual only an ARMA (p,q) process. Get the intuition on stationary (here).

fit1 <- auto.arima(DJI.ret, trace=TRUE, test="kpss",  ic="bic")

The above command gave us an ARIMA(1,0,1) or ARMA(1,1) model. Usually, the high frequency data fails the basic econometric assumption of normality, non autocorrelation and heteroskedasticity.

STEP-5: Test for ARCH Effect

However to proceed ahead with GARCH model, the residual of ARIMA must need to have the ARCH effect. To test the ARCH, lets perform the Ljung-Box test on the first 12 lags of the squared residuals of the best ARIMA model under the null hypothesis of no ARCH effects.

Box.test(fit1$residuals^2,lag=12, type="Ljung-Box")

If p-value of Ljung-Box test is smaller than 5% level of significance then there exist the ARCH effect which shows the green light to proceed ahead to GARCH.

STEP-6: Developing a GARCH Model

Instead to fit a GARCH(P,Q), in this blog, I will only fit the GARCH(1,1) for the sake of simplicity. For GARCH theory, check my blog (here).  Lets specify object called res_garch11_spec in which we want to develop a GARCH(1,1) on ARIMA(1,0,1) or in our case ARMA(1,1).

res_garch11_spec <- ugarchspec(variance.model = list(garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 1)))

Then, lets fit such specification in the data.

res_garch11_fit <- ugarchfit(spec = res_garch11_spec, data = DJI.ret)
res_garch11_fit

STEP-7: Backtesting the risk model

Once we have fitted the GARCH which modeled the risk. We would like to check the model performance by performing historical backtest. For that, we can compare the estimated VaR (value at risk) with the actual return over the period. If the return is more negative than the VaR, we have a VaR exceedance. In our case, a VaR exceedance should only occur in 1% of the cases  if we specified a 99% confidence level. VaR at 1% shows what 1% probability of my extreme loss.

For that lets implement a rolling window (refit.window = "moving") with start period of the backtest of 120 i.e n.start = 120, reestimated each day (refit.every = 1) implementing hybrid solver (solver = "hybrid") and calculate the VaR (calculate.VaR = TRUE) at the 99% VaR tail level (VaR.alpha = 0.01). And finally lets keep the coefficients (keep.coef = TRUE). There can be chances of model not to be converged so lets implement the ctrl = list(tol = 1e-7, delta = 1e-9) command prior to estimate by a rolling window.

ctrl = list(tol = 1e-7, delta = 1e-9)
res_garch11_roll <- ugarchroll(res_garch11_spec, DJI.ret, n.start = 120, refit.every = 1, refit.window = "moving", solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = 0.01, keep.coef = TRUE, solver.control = ctrl, fit.control = list(scale = 1))

STEP-8: Kupiec's Unconditional Coverage and the Christoffersen Test

We can examine the backtesting report using the report function. By specifying the type argument as VaR, the function executes the unconditional and conditional coverage tests for exceedances. VaR.alpha is the tail probability and conf.level is the confidence level on which the conditional coverage hypothesis test will be based.

report(res_garch11_roll, type = "VaR", VaR.alpha = 0.01, conf.level = 0.99)

To compare the number of expected versus actual exceedances given the tail probability of VaR, we can see the Kupiec's unconditional coverage, while the Christoffersen test tests a joint test of the unconditional coverage and the independence of the exceedances. If the actual exceedances are greater than the expected exceedances then we can't reject the null hypothesis that the exceedances are correct and independent.

STEP-9: Plotting the VaR Limits

A plot of the backtesting performance can also be generated easily using following command and selecting 2.

plot(res_garch11_fit)
selection: 2

Fig-3

Here, in the figure the returns of the data (blue) hits the 1% Var (Red) 47 times compare to 13 times expected. Hence we can state the rejection of the null hypothesis that the exceedances are correct and independent. The p-values of Kupiec's unconditional coverage is also less than 5% suggests same. Hence model shouldn't be used in real life. So forecasting is very hazardious however we will forecast for the sake of learning. To mitigate the problem, probably, we should find best GARCH(P,Q) model rather than GARCH(1,1) or may be we need to model other extensions of GARCH.

STEP-10: Forecasting Risk and VaR

For forecasting we can implement following command. Here we just forecast for 12 time period ahead.

res_garch11_fcst <- ugarchforecast(res_garch11_fit, n.ahead = 12)
res_garch11_fcst

The forecast values for 12 period are given below.

 Series  Sigma
T+1  0.09602 0.7462
T+2  0.06433 0.7537
T+3  0.06650 0.7609
T+4  0.06635 0.7677
T+5  0.06636 0.7742
T+6  0.06636 0.7804
T+7  0.06636 0.7863
T+8  0.06636 0.7919
T+9  0.06636 0.7973
T+10 0.06636 0.8025
T+11 0.06636 0.8074
T+12 0.06636 0.8121

The forecasts for the volatility are given by sigma and 1% VaR is given just by the side for each period. For the period T+3, the expected volatility is 0.0665 and the 1% VaR is 0.7609.

Reference:
I used and modified the R codes of book "Introduction to R for Quantitative Finance" which is combinely written by Gergely Daróczi, Michael Puhle, Edina Berlinger, Péter Csóka, Dániel Havran, Márton Michaletzky, Zsolt Tulassay, Kata Váradi and Agnes Vidovics-Dancs

26 comments:

  1. Thanks, Shishir, for this tutorial. Now if you have to write out the equation of the GARCH model based on the output from R, what will it be?

    ReplyDelete
    Replies
    1. http://shishirshakya.blogspot.com/2015/07/generalized-autoregressive-conditional.html

      Delete
  2. This tutorial is really helpful. But is it possible to create plot for series 5% VaR Limits?

    ReplyDelete
  3. Yes dear Ola, use the following command and select 2 to generate 5% VaR.
    plot(res_garch11_fit)
    selection: 2

    ReplyDelete
    Replies
    1. I did it as you said and "selection 2" is "Series with 1% VaR Limits". In that command there's no option to choose "Series with 5% VaR Limits".

      Delete
    2. Same problem as Ola. Can only select "Series with 1% VaR Limits", and not other option. Can you tell how to create a 5% plot?

      Delete
  4. Hey, thanks for your tutorial, it is really very helpful. But I have a question about the 'auto.arima'. when I run it the result shows that 'Best model: ARIMA(0,0,0) with non-zero mean'. I use the same data as you used. May I ask why?

    ReplyDelete
  5. Why did you set the n.Start in a STEP-7 to 120? Is there any good practice like a certain percentage of the observations? Thanks for the answer and for your great work:)

    ReplyDelete
    Replies
    1. Frankly, there can be lot of debate. If I want to generate good out of sample forecast say based upon MAPE, I would iterate and define that sweet spot for the split.

      Delete
    2. Frankly, there can be lot of debate. If I want to generate good out of sample forecast say based upon MAPE, I would iterate and define that sweet spot for the split.

      Delete
  6. I have few questions to ask.. i tried with the code for my data and i am having few errors.
    when i tried to solve the issue. I didnt get any solution for that.

    1.what if the p-value is greater in Box.test?
    2.Warning message:
    In .makefitmodel(garchmodel = "sGARCH", f = .sgarchLLH, T = T, m = m, :
    rugarch-->warning: failed to invert hessian

    Please help me to solve the issue.

    ReplyDelete
  7. Is it "Hence we can state the rejection of the null hypothesis that the exceedances are correct and independent." or "Hence we can reject the null hypothesis that the exceedances are correct and independent"? I struggled to understand that statement with the context of what is happening.And again "Here, in the figure the returns of the data (blue) hits the 1% Var (Red) 47 times compare to 13 times expected." It should be clearly stated that the values 47 and 13 are coming from the report.

    ReplyDelete
  8. Thanks for your introduction.

    After estimate the VaR, can we extract the data to Excel?
    I need to estimate the daily CVaR (expected shortfall).

    ReplyDelete
  9. I can not use the function. The following message appears: Error in get.hist.quote(instrument = "^gspc", start = "1998-01-01", quote = "Close"): cannot open URL 'https://ichart.finance.yahoo.com/table.csv?s=^gspc&a=0&b=01&c=1998&d=4&e=21&f=2017&g=d&q=q&y=0&z=^gspc&x=.csv'

    Do you can help me?

    ReplyDelete
  10. Yahoo have changes some policy that the reason the get.hist.quote of tseries package not working. I emailed to Adrian Trapletti (author of tseries) and he said Kurt Hornik (maintainer of tseries) is fixing the problem. However, you can try the Quandl package or quantmod package.

    ReplyDelete
  11. The tseries changed the "AdjClose" to "Adjusted".

    DJI.data = get.hist.quote(instrument="^DJI", start="2010-01-01", end="2015-07-23", quote="Adjusted", provider="yahoo", compression="d", retclass="zoo")

    ReplyDelete
  12. Error Correction Model (ECM) Panel Data EVIEWS 9
    https://www.youtube.com/watch?v=ZgCwrb6kI7w
    video Introduce the concept of an Error Correction Model (ECM) Panel Data EVIEWS 9.
    WhatsApp : +6285227746673
    PIN BB : D04EBECB
    IG : @olahdatasemarang

    ReplyDelete
  13. if I've run the back-test for two different garch models how can i compare them to see which is better

    ReplyDelete
  14. Hi
    In this your video you got outputs but you didn't interpret, if you can interpret it is really helpful to me.
    can you please send me the GARCH model estimation & model diagnostic link to my mail or post it

    ReplyDelete
  15. Hi.
    Is there any possibility to fit SARIMA-GARCH model in "rugarch" package?
    thank you.

    ReplyDelete
  16. Hi, it's very useful. However, I want to ask the question that why it could be a problem if using log return instead of simple return to calculate VAR as well as ES in parametric approach?

    ReplyDelete
  17. Thanks ya, artikel sangat membantu dalam menyelesaikan tugas perkuliahan tentang Generalized AutoRegressive Conditional Heteroskedastisitas (GARCH). Kunjungi juga ya MAKALAH GARCH  

    ReplyDelete
  18. Hi,Very great tutorial and explanation.I want to ask the question that ugarchforecast predicts volatility for next 12 time periods.How can we interprete this volatility.For the period T+3, the expected volatility is 0.0665 and the 1% VaR is 0.7609.How can we relate with DJI return for time periord T+3. What is the corresponding return value?How should we interprete volatility and cooresponding return?
    Thanks

    ReplyDelete
  19. if we want to backtest another method like historic simulation or conditionel extreme value theory ?we did what??

    ReplyDelete