Library and Workspace
library(readxl)
library(car)
## Loading required package: carData
library(fUnitRoots)
## Warning: package 'fUnitRoots' was built under R version 4.4.1
library(urca)
##
## Attaching package: 'urca'
## The following objects are masked from 'package:fUnitRoots':
##
## punitroot, qunitroot, unitrootTable
library(vars)
## Warning: package 'vars' was built under R version 4.4.1
## Loading required package: MASS
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.4.1
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: lmtest
library(AER)
## Warning: package 'AER' was built under R version 4.4.1
## Loading required package: survival
library(lmtest)
library(rugarch)
## Warning: package 'rugarch' was built under R version 4.4.1
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
library(rmgarch)
## Warning: package 'rmgarch' was built under R version 4.4.1
load("D:/IU/MA/3/Summer/R/Dataset/currencies.RData")
load("D:/IU/MA/3/Summer/R/Dataset/sp500.RData")
Volatility Modelling
Estimating GARCH Models
spec = ugarchspec (mean.model = list(armaOrder =c(0,0)),variance.model = list(garchOrder =c(1,1),model ="sGARCH"))
ugarchfit(spec, data = currencies$rjpy)
##
## *---------------------------------*
## * 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.006403 0.004903 1.3059 0.19158
## omega 0.001640 0.000216 7.6027 0.00000
## alpha1 0.036849 0.001690 21.8038 0.00000
## beta1 0.956418 0.001184 808.0811 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.006403 0.005947 1.0768 0.28157
## omega 0.001640 0.000608 2.6964 0.00701
## alpha1 0.036849 0.003673 10.0336 0.00000
## beta1 0.956418 0.000656 1457.4040 0.00000
##
## LogLikelihood : -4344.782
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.2180
## Bayes 1.2218
## Shibata 1.2180
## Hannan-Quinn 1.2193
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 181.0 0
## Lag[2*(p+q)+(p+q)-1][2] 184.5 0
## Lag[4*(p+q)+(p+q)-1][5] 189.1 0
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 20.79 5.114e-06
## Lag[2*(p+q)+(p+q)-1][5] 25.75 4.923e-07
## Lag[4*(p+q)+(p+q)-1][9] 31.56 1.634e-07
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 3.959 0.500 2.000 0.046614
## ARCH Lag[5] 8.886 1.440 1.667 0.012264
## ARCH Lag[7] 13.405 2.315 1.543 0.002842
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.1106
## Individual Statistics:
## mu 0.1077
## omega 0.3804
## alpha1 0.4109
## beta1 0.5065
##
## 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 1.1641 2.444e-01
## Negative Sign Bias 6.0959 1.145e-09 ***
## Positive Sign Bias 0.9241 3.555e-01
## Joint Effect 39.5060 1.356e-08 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 1911 0
## 2 30 2522 0
## 3 40 3038 0
## 4 50 3460 0
##
##
## Elapsed time : 0.6470509
EGARCH and GJR Models
spec = ugarchspec (mean.model = list(armaOrder =c(0,0)),variance.model = list(garchOrder =c(1,1),model ="eGARCH"))
ugarchfit(spec, data = currencies$rjpy)
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.003186 0.004741 0.67207 0.50154
## omega -0.011004 0.003004 -3.66269 0.00025
## alpha1 -0.028417 0.001212 -23.43837 0.00000
## beta1 0.986437 0.002677 368.53083 0.00000
## gamma1 0.102903 0.007283 14.12952 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.003186 0.006297 0.50595 0.612894
## omega -0.011004 0.009325 -1.18000 0.238001
## alpha1 -0.028417 0.009386 -3.02763 0.002465
## beta1 0.986437 0.007030 140.31135 0.000000
## gamma1 0.102903 0.026863 3.83062 0.000128
##
## LogLikelihood : -4323
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.2122
## Bayes 1.2170
## Shibata 1.2122
## Hannan-Quinn 1.2138
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 182.6 0
## Lag[2*(p+q)+(p+q)-1][2] 185.7 0
## Lag[4*(p+q)+(p+q)-1][5] 189.8 0
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 16.41 5.109e-05
## Lag[2*(p+q)+(p+q)-1][5] 22.55 3.845e-06
## Lag[4*(p+q)+(p+q)-1][9] 29.87 4.736e-07
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 4.453 0.500 2.000 0.0348425
## ARCH Lag[5] 10.646 1.440 1.667 0.0045694
## ARCH Lag[7] 16.477 2.315 1.543 0.0004957
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.8241
## Individual Statistics:
## mu 0.09764
## omega 0.70287
## alpha1 0.17568
## beta1 0.57183
## gamma1 0.68365
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.350 1.769e-01
## Negative Sign Bias 5.095 3.575e-07 ***
## Positive Sign Bias 1.124 2.610e-01
## Joint Effect 27.335 5.008e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 1749 0
## 2 30 2308 0
## 3 40 2827 0
## 4 50 3225 0
##
##
## Elapsed time : 1.079286
spec = ugarchspec (mean.model = list(armaOrder =c(0,0)),variance.model = list(garchOrder =c(1,1),model ="gjrGARCH"))
ugarchfit(spec, data = currencies$rjpy)
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.003522 0.004913 0.7169 0.473438
## omega 0.002037 0.000333 6.1231 0.000000
## alpha1 0.027579 0.003291 8.3791 0.000000
## beta1 0.950473 0.004008 237.1241 0.000000
## gamma1 0.026850 0.005349 5.0195 0.000001
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.003522 0.005938 0.59322 0.553036
## omega 0.002037 0.000764 2.66594 0.007677
## alpha1 0.027579 0.006218 4.43533 0.000009
## beta1 0.950473 0.006388 148.79911 0.000000
## gamma1 0.026850 0.010414 2.57816 0.009933
##
## LogLikelihood : -4329.648
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.2140
## Bayes 1.2188
## Shibata 1.2140
## Hannan-Quinn 1.2157
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 185.1 0
## Lag[2*(p+q)+(p+q)-1][2] 188.2 0
## Lag[4*(p+q)+(p+q)-1][5] 192.4 0
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 15.96 6.468e-05
## Lag[2*(p+q)+(p+q)-1][5] 23.44 2.163e-06
## Lag[4*(p+q)+(p+q)-1][9] 30.42 3.356e-07
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 5.856 0.500 2.000 0.0155276
## ARCH Lag[5] 12.548 1.440 1.667 0.0015551
## ARCH Lag[7] 17.377 2.315 1.543 0.0002948
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.5211
## Individual Statistics:
## mu 0.05841
## omega 0.47122
## alpha1 0.47883
## beta1 0.63269
## gamma1 0.35973
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.9567 3.388e-01
## Negative Sign Bias 5.0233 5.202e-07 ***
## Positive Sign Bias 1.3087 1.907e-01
## Joint Effect 27.4114 4.827e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 1795 0
## 2 30 2328 0
## 3 40 2892 0
## 4 50 3214 0
##
##
## Elapsed time : 1.53042
GARCH-M Estimation
meanmodel = list(armaOrder =c(0,0),archm =T, archpow =2)
varmodel = list(garchOrder =c(1,1),model ="sGARCH")
spec = ugarchspec(mean.model = meanmodel , variance.model = varmodel)
ugarchfit(spec , data = currencies$rjpy)
##
## *---------------------------------*
## * 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.011268 0.010089 1.11679 0.26409
## archm -0.027702 0.050193 -0.55192 0.58101
## omega 0.001626 0.000156 10.44456 0.00000
## alpha1 0.036768 0.000379 96.99834 0.00000
## beta1 0.956565 0.001895 504.80735 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.011268 0.012151 0.92732 0.353762
## archm -0.027702 0.053756 -0.51533 0.606322
## omega 0.001626 0.000600 2.70965 0.006735
## alpha1 0.036768 0.003722 9.87807 0.000000
## beta1 0.956565 0.001787 535.35329 0.000000
##
## LogLikelihood : -4344.631
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.2182
## Bayes 1.2230
## Shibata 1.2182
## Hannan-Quinn 1.2199
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 181.0 0
## Lag[2*(p+q)+(p+q)-1][2] 184.5 0
## Lag[4*(p+q)+(p+q)-1][5] 189.1 0
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 21.08 4.414e-06
## Lag[2*(p+q)+(p+q)-1][5] 26.01 4.178e-07
## Lag[4*(p+q)+(p+q)-1][9] 31.82 1.382e-07
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 3.933 0.500 2.000 0.047344
## ARCH Lag[5] 8.835 1.440 1.667 0.012622
## ARCH Lag[7] 13.384 2.315 1.543 0.002875
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.2634
## Individual Statistics:
## mu 0.08123
## archm 0.06069
## omega 0.38441
## alpha1 0.40883
## beta1 0.50746
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.065 2.870e-01
## Negative Sign Bias 6.068 1.358e-09 ***
## Positive Sign Bias 1.001 3.169e-01
## Joint Effect 39.673 1.250e-08 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 1799 0
## 2 30 2339 0
## 3 40 2826 0
## 4 50 3161 0
##
##
## Elapsed time : 3.205062
Forecasting from GARCH Models
spec = ugarchspec(mean.model = list(armaOrder =c(2,3)),variance.model = list(garchOrder =c(1,1),model ="sGARCH"))
fit = ugarchfit(spec, data = currencies$rjpy, out.sample = 700)
static_fc = ugarchforecast(fit ,n.ahead =1, n.roll = 699)
dynamic_fc = ugarchforecast(fit ,n.ahead = 700)
par (lwd =2,cex.axis = 2)
x_axis = currencies$Date[currencies$Date >= "2016-08-03"]
plot (x_axis , static_fc@forecast$sigmaFor, type ="l", xlab ="", ylab ="", col ="blue3")
lines (x_axis , dynamic_fc@forecast$sigmaFor ,col=" brown 1")
legend ("bottomleft", legend =c("Static", "Dynamic"),col =c("blue3","brown1"),lty= 1)
par (lwd =2,cex.axis = 2)
x_axis = currencies$Date[currencies$Date >= "2016-08-03"]
plot (x_axis , static_fc@forecast$seriesFor, type ="l", xlab ="", ylab ="", col ="blue3")
lines (x_axis , dynamic_fc@forecast$seriesFor ,col=" brown 1")
legend ("bottomleft", legend =c("Static", "Dynamic"),col =c("blue3","brown1"),lty= 1)
Estimation of Multivariate GARCH Models
uspec = ugarchspec(mean.model = list(armaOrder =c(0,0)),variance.model = list(garchOrder =c(1,1),model ="sGARCH"))
mspec = multispec(replicate(3,uspec))
cccspec = cgarchspec(mspec ,VAR = F)
mod = cgarchfit(cccspec ,data = currencies[c("reur","rgbp","rjpy")])
mod
##
## *-------------------------------------------------*
## * Copula GARCH Fit *
## *-------------------------------------------------*
##
## Distribution : mvnorm
## No. of Parameters : 12
## [VAR GARCH CC] : [0+12+0]
## No. of Series : 3
## No. of Observations : 7141
## Log-Likelihood : -9464.444
## Av.Log-Likelihood : -1.325
##
## Optimal Parameters
## ---------------------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## [reur].mu -0.003892 0.004654 -0.83629 0.402993
## [reur].omega 0.000272 0.000124 2.19087 0.028461
## [reur].alpha1 0.020849 0.000967 21.56751 0.000000
## [reur].beta1 0.978108 0.000040 24503.86687 0.000000
## [rgbp].mu -0.001328 0.004365 -0.30422 0.760964
## [rgbp].omega 0.000827 0.000275 3.00633 0.002644
## [rgbp].alpha1 0.034702 0.003231 10.73921 0.000000
## [rgbp].beta1 0.961216 0.000506 1899.29985 0.000000
## [rjpy].mu 0.006403 0.005269 1.21518 0.224298
## [rjpy].omega 0.001640 0.000592 2.76890 0.005625
## [rjpy].alpha1 0.036849 0.003524 10.45584 0.000000
## [rjpy].beta1 0.956418 0.000614 1558.89610 0.000000
##
## Information Criteria
## ---------------------
##
## Akaike 2.6541
## Bayes 2.6656
## Shibata 2.6541
## Hannan-Quinn 2.6581
##
##
## Elapsed time : 7.703608
mod@mfit$Rt
## [,1] [,2] [,3]
## [1,] 1.0000000 0.6374866 0.3153877
## [2,] 0.6374866 1.0000000 0.2195692
## [3,] 0.3153877 0.2195692 1.0000000
Simulation method
Deriving Critical Values for a Dickey{Fuller Test Using
Simulation
set.seed(123456)
tnone = tconst = ttrend = NULL
for (i in 1:500){
y = c(0, cumsum(rnorm(1199)))
dy = c(NA , diff ( y ) )
lagy = c(NA , y[1:1199])
# Model with no intercept (no drift)
lmnone = lm(dy[-(1:200)] ~ 0 + lagy[-(1:200)])
tnone = c(tnone , coef(summary(lmnone))[1]/coef(summary(lmnone))[2])
# Model with intercept, no trend
lmconst = lm(dy[-(1:200)] ~ 1 + lagy[-(1:200)])
tconst = c(tconst, coef(summary(lmconst))[2,1]/coef(summary(lmconst))[2,2])
# Model with intercept and trend
lmtrend = lm(dy[-(1:200)] ~ 1 + lagy[-(1:200)] + c(1:1000))
ttrend = c(ttrend, coef(summary(lmtrend))[2,1]/ coef(summary(lmtrend))[2,2])
}
quantile(tnone, c(0.01,0.05,0.1))
## 1% 5% 10%
## -2.401095 -1.921597 -1.675099
quantile(tconst, c(0.01,0.05,0.1))
## 1% 5% 10%
## -3.651240 -2.870267 -2.571466
quantile(ttrend, c(0.01,0.05,0.1))
## 1% 5% 10%
## -3.987071 -3.366447 -3.107267
Pricing Asian Options
set.seed(123456)
# Initial values and derived constants
K = 6500; iv = 0.2652 # Simulation one
# K = 5500; iv = 0.3433 # Simulation two
s0 = 6289.7; rf = 0.0624; dy = 0.0242; ttm = 0.5; obs = 125
dt = ttm /obs
drift = (rf -dy -iv^2/2)*dt
vsqrdt = iv*dt^0.5
putval = callval = NULL
for (i in 1:25000){
random = rnorm (obs ) # create cumulative sum of random numbers
# Spot price evolution for positive and negative random numbers
spot = s0* exp( drift *c(1: obs)+ vsqrdt * cumsum ( random ))
spot_neg = s0*exp( drift *c(1:obs)+ vsqrdt * cumsum (- random ))
# Compute call values
callval = c( callval , max ( mean ( spot )-K,0)*exp (-rf*ttm ))
callval = c( callval , max ( mean ( spot_neg )-K,0)*exp (-rf*ttm))
# Compute Put values
putval = c(putval , max (K- mean ( spot ),0)*exp (-rf*ttm ))
putval = c(putval , max (K- mean ( spot_neg ),0)*exp (-rf*ttm ))
}
mean(callval)
## [1] 204.6794
mean(putval)
## [1] 349.6241