# Set all the working directory 
library(readxl)
library(dplyr)
library(xts)
library(zoo)
library(tseries)
library(forecast)
library(astsa)
library(lmtest)

#Cumulative volatility and VaR

library(quantmod)
## Warning: package 'quantmod' was built under R version 4.0.3
## Loading required package: TTR
MXindex<-getSymbols("^MXX", from="2018-01-01", to="2021-04-28",auto.assign=FALSE)[,6]
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
MXreturn<- na.omit(diff(log(MXindex)))
names(MXreturn)<-c("return")

hist(MXreturn, breaks=40)

Now I run a GARCH(1,1) model using the daily returns

library(rugarch)
## Warning: package 'rugarch' was built under R version 4.0.5
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
Especgarch1=ugarchspec(variance.model= list(model= "sGARCH", garchOrder= c(1, 1), 
submodel= NULL,  variance.targeting= FALSE), 
mean.model= list(armaOrder= c(0, 0), include.mean= TRUE, archm= FALSE, 
archpow= 0, arfima= FALSE,  external.regressors= NULL, archex= FALSE), 
distribution.model= "norm", start.pars= list(), fixed.pars= list())

#Estimación del GARCH(1,1)

garch1<- ugarchfit(spec=Especgarch1, data=MXreturn$return)
garch1
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000067    0.000323  0.20726 0.835809
## omega   0.000006    0.000002  3.13733 0.001705
## alpha1  0.143639    0.021486  6.68538 0.000000
## beta1   0.810415    0.019654 41.23353 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000067    0.000329  0.20332 0.838887
## omega   0.000006    0.000003  2.03500 0.041851
## alpha1  0.143639    0.019550  7.34725 0.000000
## beta1   0.810415    0.033569 24.14150 0.000000
## 
## LogLikelihood : 2623.694 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.3049
## Bayes        -6.2822
## Shibata      -6.3050
## Hannan-Quinn -6.2962
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic  p-value
## Lag[1]                      7.359 0.006672
## Lag[2*(p+q)+(p+q)-1][2]     7.578 0.008369
## Lag[4*(p+q)+(p+q)-1][5]     8.264 0.025327
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.5599 0.45430
## Lag[2*(p+q)+(p+q)-1][5]    3.5506 0.31564
## Lag[4*(p+q)+(p+q)-1][9]   10.2552 0.04422
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale  P-Value
## ARCH Lag[3]     3.645 0.500 2.000 0.056230
## ARCH Lag[5]     5.272 1.440 1.667 0.088905
## ARCH Lag[7]    11.592 2.315 1.543 0.007767
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  6.369
## Individual Statistics:             
## mu     0.2703
## omega  0.5169
## alpha1 0.2165
## beta1  0.3623
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.4834 0.6289    
## Negative Sign Bias  0.1837 0.8543    
## Positive Sign Bias  0.3662 0.7143    
## Joint Effect        1.4094 0.7033    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     26.83      0.10858
## 2    30     36.91      0.14864
## 3    40     56.36      0.03549
## 4    50     62.20      0.09755
## 
## 
## Elapsed time : 0.1860139
##

Now I do a prediction for the next 10 days

garch1.prediction <- ugarchforecast(garch1,n.ahead= 10)
garch1.prediction@forecast$seriesFor
##        2021-04-27
## T+1  6.690189e-05
## T+2  6.690189e-05
## T+3  6.690189e-05
## T+4  6.690189e-05
## T+5  6.690189e-05
## T+6  6.690189e-05
## T+7  6.690189e-05
## T+8  6.690189e-05
## T+9  6.690189e-05
## T+10 6.690189e-05
#I get the forecast for the daily volatility for the next 10 days
pred_vol= garch1.prediction@forecast$sigmaFor
pred_vol
##       2021-04-27
## T+1  0.008016813
## T+2  0.008204170
## T+3  0.008379016
## T+4  0.008542493
## T+5  0.008695594
## T+6  0.008839190
## T+7  0.008974046
## T+8  0.009100844
## T+9  0.009220191
## T+10 0.009332631

I can get the acumulative volatily for the 10 days

I first get the variances from the forecasted standar deviations (volatilities)

pred_variance = pred_vol^2
pred_variance
##        2021-04-27
## T+1  6.426929e-05
## T+2  6.730841e-05
## T+3  7.020791e-05
## T+4  7.297418e-05
## T+5  7.561335e-05
## T+6  7.813127e-05
## T+7  8.053350e-05
## T+8  8.282536e-05
## T+9  8.501192e-05
## T+10 8.709801e-05

Now I can ADD these 10 variances

pred_variance10= sum(pred_variance)
pred_variance10
## [1] 0.0007639732

Now that I have the 10-day variance, I can easily get the 10-day volatility:

pred_volatility10=(sqrt(pred_variance10))
pred_volatility10
## [1] 0.02764007

Now I can get the 10-day 1%VaR

p%VaR = mean(return) + t * volatility

library(moments)
## Warning: package 'moments' was built under R version 4.0.3
meanreturn = mean(MXreturn$return)
k = kurtosis(MXreturn$return)
df = 6/k+4
df
##   return 
## 4.891423
# I calculate the t value for the VaR according to 1% probability
t = qt(0.01, df)
t
## [1] -3.396388
z = qnorm(0.01)
z
## [1] -2.326348
# Now I can calculate the cumulative  1%VaR for the next 10 days: 

VaR1 <- meanreturn + t *pred_volatility10
VaR1
## [1] -0.09390613

THE -0.09390613 MEANS THAT I CAN LOOSE -9.390613% OR MORE WITH A PROBABILITY OF 1%

1 Market GARCH model

I download daily data for the stock

SP1<-getSymbols("AMXL.MX", from="2018-01-01", to="2021-04-28",auto.assign=FALSE)[,6]
sreturn<- na.omit(diff(log(SP1)))
names(sreturn)<-c("return")

hist(sreturn, breaks=40)

Now I run the model including the mean model and the GARCH(1,1) , volatility model:

Especgarch2=ugarchspec(variance.model= list(model= "sGARCH", garchOrder= c(1, 1), 
submodel= NULL,  variance.targeting= FALSE), 
mean.model= list(armaOrder= c(1, 1), include.mean= TRUE, archm= FALSE, 
archpow= 0, arfima= FALSE,  external.regressors= NULL, archex= FALSE), 
distribution.model= "norm", start.pars= list(), fixed.pars= list())

#Estimación del GARCH(1,1)

garch2<- ugarchfit(spec=Especgarch2, data=sreturn$return)
garch2
## 
## *---------------------------------*
## *          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      0.000033    0.000532  0.061914 0.950631
## ar1    -0.498611    0.320079 -1.557775 0.119287
## ma1     0.539548    0.309876  1.741176 0.081653
## omega   0.000005    0.000003  1.658450 0.097227
## alpha1  0.049666    0.008385  5.923487 0.000000
## beta1   0.931219    0.009896 94.096159 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000033    0.000479  0.068744 0.945194
## ar1    -0.498611    0.175364 -2.843285 0.004465
## ma1     0.539548    0.176364  3.059292 0.002219
## omega   0.000005    0.000007  0.718407 0.472507
## alpha1  0.049666    0.015111  3.286713 0.001014
## beta1   0.931219    0.021670 42.972248 0.000000
## 
## LogLikelihood : 2261.316 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.4280
## Bayes        -5.3939
## Shibata      -5.4281
## Hannan-Quinn -5.4149
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2249  0.6353
## Lag[2*(p+q)+(p+q)-1][5]    2.9437  0.5076
## Lag[4*(p+q)+(p+q)-1][9]    6.9205  0.1358
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  1.040e-07  0.9997
## Lag[2*(p+q)+(p+q)-1][5] 1.750e+00  0.6785
## Lag[4*(p+q)+(p+q)-1][9] 3.092e+00  0.7437
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     2.749 0.500 2.000 0.09732
## ARCH Lag[5]     2.964 1.440 1.667 0.29501
## ARCH Lag[7]     3.765 2.315 1.543 0.38168
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.3438
## Individual Statistics:              
## mu     0.10644
## ar1    0.27092
## ma1    0.23713
## omega  0.10207
## alpha1 0.10418
## beta1  0.07832
## 
## 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           0.5383 0.5905    
## Negative Sign Bias  0.3828 0.7020    
## Positive Sign Bias  0.8400 0.4012    
## Joint Effect        0.8767 0.8310    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     37.18     0.007527
## 2    30     49.61     0.009937
## 3    40     65.03     0.005547
## 4    50     81.09     0.002667
## 
## 
## Elapsed time : 0.3430421

AR1 (-0.498611) IS NEGATIVE AND SIGNIFICANT BECAUSE WHEN A STOCK GAINS A LOT IN 1 DAY, THE NEXT DAY IT GETS A NEGATIVE VALUE.ALSO MEANS THAT ONLY -49.86% OF YESTERDAY’S INFORMATION IS PASSED TO TODAY’S RETURNS. MA1 (0.539548) IS POSITIVE AND STATISTICALLY SIGNIFICANT, THIS MEANS THAT THE STOCKS OF YESTERDAY ARE CORRELATED WITH TODAY’S

BETA1 (0.931219) IS PORITIVE AND STATISTICALLY SIGNIFICANT, MEANING THAT ARROUND 93.1219% OF YESTERDAY´S RETURNS VARIANCE IS PASSED TO TODA’S RETURNS VARIANCE. ALSO AS THE VALUE IS GREATER THAN 0.9, THEN THIS MEANS THERE’S CLUSTERING OF VOLATILITY.

2 Practicing with ARCH models and Value at Risk

Using the same model you did in the previous part, do the following:

Using volatility estimations from this model, you have to calculate 5%VaR, 5%ES for each of the following 10 days in the future assuming fat-tailed distri- bution. Also, you have to estimate the cumulative 5%VaR and the cumulative 5%ES for the following 5 days.

garch2.prediction <- ugarchforecast(garch2,n.ahead= 10)

pred_vol2 = garch2.prediction@forecast$sigmaFor
pred_vol2
##      2021-04-27
## T+1  0.01246630
## T+2  0.01255536
## T+3  0.01264212
## T+4  0.01272664
## T+5  0.01280900
## T+6  0.01288927
## T+7  0.01296753
## T+8  0.01304384
## T+9  0.01311826
## T+10 0.01319084
sp2 <-getSymbols("AMXL.MX", from="2018-01-01", to="2021-04-28",auto.assign=FALSE)[,6]
AMXLreturn <- na.omit(diff(log(sp2)))
names(AMXLreturn)<-c("return")

prices<-merge(Ad(sp2),Ad(MXindex))
r <- na.omit(diff(log(prices)))
kurtosis2 = as.numeric(kurtosis(sp2))
kurtosis2
## [1] 2.643347
mean_sp = mean(AMXLreturn)
df = (6/kurtosis2) + 4
df
## [1] 6.269849
t1 = qt(0.05,df)

VaR_5 = mean_sp + t1 * pred_vol2
VaR_5
##       2021-04-27
## T+1  -0.02412866
## T+2  -0.02430041
## T+3  -0.02446769
## T+4  -0.02463067
## T+5  -0.02478949
## T+6  -0.02494428
## T+7  -0.02509519
## T+8  -0.02524233
## T+9  -0.02538583
## T+10 -0.02552580

2.1 5%Expected Shortfall

t2 = qt(0.05/2,df)
ES5 = mean_sp + t2*pred_vol2
ES5
##       2021-04-27
## T+1  -0.03027834
## T+2  -0.03049402
## T+3  -0.03070410
## T+4  -0.03090877
## T+5  -0.03110822
## T+6  -0.03130262
## T+7  -0.03149213
## T+8  -0.03167691
## T+9  -0.03185712
## T+10 -0.03203290
#cumulative 5%VaR
garch2.prediction <- ugarchforecast(garch2,n.ahead= 5)
pred_vol3 = garch2.prediction@forecast$sigmaFor
pred_vol3
##     2021-04-27
## T+1 0.01246630
## T+2 0.01255536
## T+3 0.01264212
## T+4 0.01272664
## T+5 0.01280900
#I elevate to get the standard deviations
pred_variance2 = pred_vol3^2
pred_variance2
##       2021-04-27
## T+1 0.0001554087
## T+2 0.0001576372
## T+3 0.0001598231
## T+4 0.0001619673
## T+5 0.0001640704
#Sum the standard deviations to get the cumulative 
variance = sum(pred_variance2)
variance
## [1] 0.0007989066
#Get the squared root of the sum to reach the volatility
vol = sqrt(variance)
vol
## [1] 0.02826494
#5%VaRp%VaR=mean(return)+t*volatility
Acum_VaR = mean_sp + t1 * vol
Acum_VaR
## [1] -0.05459325

At a 5% PROBABILITY THE MINIMUM LOSS THAT AMX WOULD HAVE IS ABOUT -05.459326%

#acum5% expected shortfall
Acum.es5 = mean_sp + t2*vol
Acum.es5
## [1] -0.06853646

THE EXPECTED SHORTFALL GIVE US -0.06853646 WHICH MEANS THAT THE EXOECTED AVERAGE LOSS OF AN INVESTMENT IN AMX IN ONE DAY ONCE THE PROBABILITY 5%P OF A BAD SCENARIO HAPPENS Expected Shortfall with a probability p refers to the expected average loss of an investment in ONE period (in this case, in one day) once the probability p of a bad scenario happens.