# http://www.r-bloggers.com/high-frequency-garch-the-multiplicative-component-garch-mcsgarch-model/
library(rugarch)
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
Sys.setenv(TZ = 'GMT')
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
R_i= read.csv("C:/Users/Asadel/Downloads/GARCH MODEL REGRESI/C_2008_1minret.csv")
R_i = xts(R_i[, 2], as.POSIXct(R_i[, 1]))
getSymbols('C', from = '2000-01-01')
## As of 0.4-0, 'getSymbols' uses env=parent.frame() and
## auto.assign=TRUE by default.
##
## This behavior will be phased out in 0.5-0 when the call will
## default to use auto.assign=FALSE. getOption("getSymbols.env") and
## getOptions("getSymbols.auto.assign") are now 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 more details.
## [1] "C"
C = adjustOHLC(C, use.Adjusted = TRUE)
R_d = ROC(Cl(C), na.pad = FALSE)
par(cex.main = 0.85, col.main = 'black')
acf(abs(as.numeric(R_i)), lag.max = 4000, main = '1-min absolute returns\nCitigroup (2008 Jan-Feb)', cex.lab = 0.8)
###
# Find the unique days in the intraday sample
n = length(unique(format(index(R_i), '%Y-%m-%d')))
# define a daily spec
spec_d = ugarchspec(mean.model = list(armaOrder = c(1, 1)), variance.model = list(model = 'eGARCH', garchOrder = c(2, 1)), distribution = 'nig')
# use the ugarchroll method to create a rolling forecast for the period in
# question:
roll = ugarchroll(spec_d, data = R_d['/2008-02-29'], forecast.length = n, refit.every = 5, refit.window = 'moving', moving.size = 2000, calculate.VaR = FALSE)
# extract the sigma forecast
df = as.data.frame(roll)
f_sigma = as.xts(df[, 'Sigma', drop = FALSE])
# now estimate the intraday model
spec = ugarchspec(mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), variance.model = list(model = 'mcsGARCH'), distribution = 'nig')
# DailyVar is the required xts object of the forecast daily variance
fit = ugarchfit(data = R_i, spec = spec, DailyVar = f_sigma^2)
###
ep <- axTicksByTime(fit@model$DiurnalVar)
par(mfrow = c(4, 1), mar = c(2.5, 2.5, 2, 1))
plot(as.numeric(fit@model$DiurnalVar^0.5), type = 'l', main = 'Sigma[Diurnal]', col = 'tomato1', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(as.numeric(fit@model$DailyVar^0.5), type = 'l', main = 'Sigma[Daily-Forecast]', col = 'tomato2', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(fit@fit$q, type = 'l', main = 'Sigma[Stochastic]', col = 'tomato3', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(as.numeric(sigma(fit)), type = 'l', main = 'Sigma[Total]', col = 'tomato4', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
##
# create the interval
interval = format(seq(as.POSIXct('2008-01-02 09:31:00'), as.POSIXct('2008-01-02 16:00:00'), by = 'min'), '%H:%M:%S')
#
ForcTime = ftseq(T0 = as.POSIXct('2008-02-29 16:00:00'), length.out = 390 * 2 + 1, by = 'mins', interval = interval)
tail(ForcTime, 25)
## [1] "2008-03-04 15:37:00 GMT" "2008-03-04 15:38:00 GMT"
## [3] "2008-03-04 15:39:00 GMT" "2008-03-04 15:40:00 GMT"
## [5] "2008-03-04 15:41:00 GMT" "2008-03-04 15:42:00 GMT"
## [7] "2008-03-04 15:43:00 GMT" "2008-03-04 15:44:00 GMT"
## [9] "2008-03-04 15:45:00 GMT" "2008-03-04 15:46:00 GMT"
## [11] "2008-03-04 15:47:00 GMT" "2008-03-04 15:48:00 GMT"
## [13] "2008-03-04 15:49:00 GMT" "2008-03-04 15:50:00 GMT"
## [15] "2008-03-04 15:51:00 GMT" "2008-03-04 15:52:00 GMT"
## [17] "2008-03-04 15:53:00 GMT" "2008-03-04 15:54:00 GMT"
## [19] "2008-03-04 15:55:00 GMT" "2008-03-04 15:56:00 GMT"
## [21] "2008-03-04 15:57:00 GMT" "2008-03-04 15:58:00 GMT"
## [23] "2008-03-04 15:59:00 GMT" "2008-03-04 16:00:00 GMT"
## [25] "2008-03-05 09:31:00 GMT"
## [1] '2008-03-04 15:37:00 GMT' '2008-03-04 15:38:00 GMT'
## [3] '2008-03-04 15:39:00 GMT' '2008-03-04 15:40:00 GMT'
## [5] '2008-03-04 15:41:00 GMT' '2008-03-04 15:42:00 GMT'
## [7] '2008-03-04 15:43:00 GMT' '2008-03-04 15:44:00 GMT'
## [9] '2008-03-04 15:45:00 GMT' '2008-03-04 15:46:00 GMT'
## [11] '2008-03-04 15:47:00 GMT' '2008-03-04 15:48:00 GMT'
## [13] '2008-03-04 15:49:00 GMT' '2008-03-04 15:50:00 GMT'
## [15] '2008-03-04 15:51:00 GMT' '2008-03-04 15:52:00 GMT'
## [17] '2008-03-04 15:53:00 GMT' '2008-03-04 15:54:00 GMT'
## [19] '2008-03-04 15:55:00 GMT' '2008-03-04 15:56:00 GMT'
## [21] '2008-03-04 15:57:00 GMT' '2008-03-04 15:58:00 GMT'
## [23] '2008-03-04 15:59:00 GMT' '2008-03-04 16:00:00 GMT'
## [25] '2008-03-05 09:31:00 GMT'
##
fit2 = ugarchfit(data = R_i, spec = spec, DailyVar = f_sigma^2, out.sample = 300)
# won't supply DailyVar to get an error
#forc = ugarchforecast(fit2, n.ahead = 10, n.roll = 299)
##
fit_d = ugarchfit(spec_d, data = R_d['2002/2008-02-29'])
forc_d = ugarchforecast(fit_d, n.ahead = 1)
f_sigma = xts(as.numeric(sigma(forc_d)), as.POSIXct('2008-03-03'))
# intraday forecast
forc = ugarchforecast(fit2, n.ahead = 10, n.roll = 299, DailyVar = f_sigma^2)
show(forc)
##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: mcsGARCH
## Horizon: 10
## Roll Steps: 299
## Out of Sample: 10
##
## 0-roll forecast [T0=2008-02-29 11:00:00]:
## Series Sigma[Total] Sigma[Stochastic]
## T+1 -7.414e-06 0.0015168 0.8599
## T+2 -7.399e-06 0.0056955 0.8616
## T+3 -7.393e-06 0.0055482 0.8632
## T+4 -7.391e-06 0.0058752 0.8649
## T+5 -7.390e-06 0.0063205 0.8665
## T+6 -7.389e-06 0.0013059 0.8681
## T+7 -7.389e-06 0.0012863 0.8698
## T+8 -7.389e-06 0.0011245 0.8713
## T+9 -7.389e-06 0.0008202 0.8729
## T+10 -7.389e-06 0.0009282 0.8745
##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: mcsGARCH
## Horizon: 10
## Roll Steps: 299
## Out of Sample: 10
##
## 0-roll forecast [T0=2008-02-29 11:00:00]:
## Series Sigma[Total] Sigma[Stochastic]
## T+1 -7.681e-06 0.0015132 0.8702
## T+2 -7.664e-06 0.0057046 0.8718
## T+3 -7.657e-06 0.0055551 0.8734
## T+4 -7.654e-06 0.0058834 0.8750
## T+5 -7.653e-06 0.0063295 0.8766
## T+6 -7.653e-06 0.0013036 0.8781
## T+7 -7.653e-06 0.0012846 0.8797
## T+8 -7.653e-06 0.0011227 0.8812
## T+9 -7.652e-06 0.0008177 0.8827
## T+10 -7.652e-06 0.0009259 0.8842
#############################################################OK#########################################
T0 = tail(index(R_i), 1)
# model$dtime contains the set of unique interval points in the dataset
# (and available from all rugarch objects for the mcsGARCH model)
# model$dvalues contains the diurnal component for each interval
ftime = ftseq(T0, length.out = 10000, by = fit@model$modeldata$period, interval = fit@model$dtime)
dtime = unique(format(ftime, '%Y-%m-%d'))
#
sim_d = ugarchsim(fit_d, n.sim = length(dtime), m.sim = 1)
var_sim = xts(as.matrix(sigma(sim_d)^2), as.POSIXct(dtime))
sim = ugarchsim(fit, n.sim = 10000, n.start = 0, m.sim = 1, DailyVar = var_sim, rseed = 10)
#
ep < - axTicksByTime(sim@simulation$DiurnalVar)
## Warning in ep < -axTicksByTime(sim@simulation$DiurnalVar): longer object
## length is not a multiple of shorter object length
## Mar 03 09:31 Mar 04 09:31 Mar 05 09:31 Mar 06 09:31 Mar 07 09:31
## FALSE FALSE FALSE FALSE FALSE
## Mar 10 09:31 Mar 11 09:31 Mar 12 09:31 Mar 13 09:31 Mar 14 09:31
## FALSE FALSE FALSE FALSE FALSE
## Mar 17 09:31 Mar 18 09:31 Mar 19 09:31 Mar 20 09:31 Mar 21 09:31
## FALSE FALSE FALSE FALSE FALSE
## Mar 24 09:31 Mar 25 09:31 Mar 26 09:31 Mar 27 09:31 Mar 28 09:31
## FALSE FALSE FALSE FALSE FALSE
## Mar 31 09:31 Apr 01 09:31 Apr 02 09:31 Apr 03 09:31 Apr 04 09:31
## FALSE FALSE FALSE FALSE FALSE
## Apr 07 09:31 Apr 07 13:40
## FALSE FALSE
par(mfrow = c(4, 1), mar = c(2.5, 2.5, 2, 1))
plot(as.numeric(sim@simulation$DiurnalVar^0.5), type = 'l', main = 'Sigma[Diurnal]', col = 'tomato1', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(as.numeric(sim@simulation$DailyVar^0.5), type = 'l', main = 'Sigma[Daily-Simulated]', col = 'tomato2', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(sim@simulation$qSim[, 1], type = 'l', main = 'Sigma[Stochastic]', col = 'tomato3', xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
plot(as.numeric(sigma(sim)), type = 'l', main = 'Sigma[Total]', col = 'tomato4',
xaxt = 'n', ylab = 'sigma', xlab = ' ')
axis(1, at = ep, labels = names(ep), tick = TRUE)
grid()
##
n = length(index(R_d['2008-01-01/2008-03-01']))
spec_d = ugarchspec(mean.model = list(armaOrder = c(1, 1)), variance.model = list(model = 'sGARCH'), distribution = 'std')
roll = ugarchroll(spec_d, data = R_d['/2008-02-29'], forecast.length = n, refit.every = 5, refit.window = 'moving', moving.size = 2000, calculate.VaR = FALSE)
df = as.data.frame(roll)
f_sigma = as.xts(df[, 'Sigma', drop = FALSE])
spec = ugarchspec(mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), variance.model = list(model = 'mcsGARCH'), distribution = 'std')
roll = ugarchroll(spec, data = R_i, DailyVar = f_sigma^2, forecast.length = 3000, refit.every = 390, refit.window = 'moving', moving.size = 3000, calculate.VaR = TRUE)
# Generate the 1% VaR report
report(roll)
## VaR Backtest Report
## ===========================================
## Model: mcsGARCH-std
## Backtest Length: 3000
## Data:
##
## ==========================================
## alpha: 1%
## Expected Exceed: 30
## Actual VaR Exceed: 33
## Actual %: 1.1%
##
## Unconditional Coverage (Kupiec)
## Null-Hypothesis: Correct Exceedances
## LR.uc Statistic: 0.294
## LR.uc Critical: 3.841
## LR.uc p-value: 0.588
## Reject Null: NO
##
## Conditional Coverage (Christoffersen)
## Null-Hypothesis: Correct Exceedances and
## Independence of Failures
## LR.cc Statistic: NaN
## LR.cc Critical: 5.991
## LR.cc p-value: NaN
## Reject Null: NA
## VaR Backtest Report
## ===========================================
## Model: mcsGARCH-std
## Backtest Length: 3000
## ==========================================
## alpha: 1%
## Expected Exceed: 30
## Actual VaR Exceed: 33
## Actual %: 1.1%
##
## Unconditional Coverage (Kupiec)
## Null-Hypothesis: Correct Exceedances
## LR.uc Statistic: 0.294
## LR.uc Critical: 3.841
## LR.uc p-value: 0.588
## Reject Null: NO
##
## Conditional Coverage (Christoffersen)
## Null-Hypothesis: Correct Exceedances and
## Independence of Failures
## LR.cc Statistic: 1.028
## LR.cc Critical: 5.991
## LR.cc p-value: 0.598
## Reject Null: NO
##
D = as.POSIXct(rownames(roll@forecast$VaR))
VaRplot(0.01, actual = xts(roll@forecast$VaR[, 3], D), VaR = xts(roll@forecast$VaR[,1], D))
