################################################################# ## "COGARCH(p,q): Simulation and Inference with Yuima package" ## ################################################################# library("yuima") ################################################### ## Section 3. Simulation of a COGARCH(p,q) model ## ################################################### ## Chunk 1: Model Definition ## model1 <- setCogarch(p = 1, q = 1, work = FALSE, measure = list(df = "rvgamma(z, 1, sqrt(2), 0, 0)"), measure.type = "code", Cogarch.var = "G", V.var = "v", Latent.var = "x", XinExpr = TRUE) ## Chunk 2: Values of model parameters ## param1 <- list(a1 = 0.038, b1 = 301, a0 = 0.01, x01 = 0) ## Chunk 3: Simulated sample path, Fig. 1 ## Terminal1 <- 5 n1 <- 750 samp1 <- setSampling(Terminal = Terminal1, n = n1) set.seed(123) sim1 <- simulate(model1, sampling = samp1, true.parameter = param1, method = "euler") # Code for generation Figure 1 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(sim1, main = "Sample Path of a VG COGARCH(1,1) model with Euler scheme") ## Chunk 4: Simulated sample path, Fig. 2 ## set.seed(123) sim2 <- simulate(model1, sampling = samp1, true.parameter = param1, method = "mixed") # Code for generation Figure 2 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(sim2, main = "Sample Path of a VG COGARCH(1,1) model with mixed scheme") ################################## ## Section 6. Numerical results ## ################################## ##################################################################### ## Subsect. 6.1. Simulation and Estimation of a COGARCH(p,q) model ## ##################################################################### ## Chunk 1: values of model parameters CP-COGARCH(1,1)## numSeed <- 200 param.cp <- list(a1 = 0.038, b1 = 0.053, a0 = 0.04 / 0.053, x01 = 50.33) ## Chunk 2: model definition CP-COGARCH(1,1)## mod.cp <- setCogarch(p = 1, q = 1, work = FALSE, measure = list(intensity = "1", df = list("dnorm(z, 0, 1)")), measure.type = "CP", Cogarch.var = "g", V.var = "v", Latent.var = "x", XinExpr = TRUE) ## Chunk 3: Sample path of a CP COGARCH(1,1) model ## Term <- 1600 num <- 24000 set.seed(numSeed) samp.cp <- setSampling(Terminal = Term, n = num) sim.cp <- simulate(mod.cp, true.parameter = param.cp, sampling = samp.cp, method = "euler") ## Chunk 4: Fig. 4 ## par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(sim.cp, main = "simulated COGARCH(1,1) model driven by a Compound Poisson process") ## Chunk 5: Estimation of a CP-COGARCH(1,1) model ## res.cp <- gmm(sim.cp, start = param.cp, objFun = "L2", Est.Incr = "Incr") ## Chunk 6: Summary 1 ## summary(res.cp) ## Chunk 7: Fig.5 ## par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(res.cp, main = "CP Increments of a COGARCH(1,1) model") ## Chunk 8: Trajectory reconstruction through estimated increments. CP-COGARCH(1,1) ## traj.cp <- simulate(res.cp) par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(traj.cp, main = "estimated COGARCH(1,1) driven by compound poisson process") ## Chunk 9: Model Definition VG-COGARCH(1,1) ## param.VG <- list(a1 = 0.038, b1 = 0.053, a0 = 0.04 / 0.053, x01 = 50.33) cog.VG <- setCogarch(p = 1, q = 1, work = FALSE, measure = list(df = "rvgamma(z, 1, sqrt(2), 0, 0)"), measure.type = "code", Cogarch.var = "y", V.var = "v", Latent.var = "x", XinExpr = TRUE) ## Chunk 10: Sample Path VG-COGARCH(1,1) model Fig. 7 ## set.seed(numSeed) samp.VG <- setSampling(Terminal = Term, n = num) sim.VG <- simulate(cog.VG, true.parameter = param.VG, sampling = samp.VG, method = "euler") # Code for generation Figure 7 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(sim.VG, main = "simulated COGARCH(1,1) model driven by a Variance Gamma process") ## Chunk 11: Estimation VG-COGARCH(1,1) Fig. 8 ## res.VG <- gmm(sim.VG, start = param.VG, Est.Incr = "Incr") summary(res.VG) # Code for generation Figure 8 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(res.VG, main = "VG Increments of a COGARCH(1,1) model") ## Chunk 12: Reconstruction of VG-COGARCH(1,1) trajectory with estimated increments Fig. 9 ## traj.VG <- simulate(res.VG) # Code for generation Figure 9 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(traj.VG, main = "estimated COGARCH(1,1) model driven by Variance Gamma process") ## Chunk 13: Definition CP-COGARCH(2,1) model ## param.cp2 <- list(a0 = 0.5, a1 = 0.1, b1 =1.5, b2 = 0.5, x01 = 2.5, x02 = 0) mod.cp2 <- setCogarch(p = 1, q = 2, work = FALSE, measure = list(intensity = "1",df = list("dnorm(z,0,1)")), measure.type = "CP", Cogarch.var = "y", V.var = "v", Latent.var = "x", XinExpr = TRUE) ## Chunk 14: Simulation of CP-COGARCH(2,1) model Fig. 10 ## samp.cp2 <- setSampling(Terminal = Term, n = num) set.seed(numSeed) sim.cp2 <- simulate(mod.cp2, sampling = samp.cp2, true.parameter = param.cp2, method="euler") # Code for generation Figure 10 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(sim.cp2, main = "simulated COGARCH(2,1) model driven by a Compound Poisson process") ## Chunk 15: Estimation of CP-COGARCH(2,1) Fig. 11 ## res.cp2 <- gmm(yuima = sim.cp2, start = param.cp2, Est.Incr = "IncrPar") summary(res.cp2) # Code for generation Figure 11 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(res.cp2, main = "Compound Poisson Increment of a COGARCH(2,1) model") ## Chunk 16: Reconstruction of CP-COGARCH(2,1) sample path with ## estimated increments Fig. 12 ## traj.cp2 <- simulate(res.cp2) # Code for generation Figure 12 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(traj.cp2, main = "estimated COGARCH(2,1) model driven by a Compound Poisson process") ############################# ## Subsect. 6.2. Real Data ## ############################# ## Chunk 1: Market Data Fig. 13 ## library("yuima") data("LogSPX", package = "yuima") allObs <- Data$allObs obsinday <- Data$obsinday logdayprice <- Data$logdayprice dimData <- length(allObs) # Code for generation Figure 13 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 3,2)) plot(x = 1:dimData / obsinday, y = allObs, main = "Intraday 5 min Log SPX Index", type = "l", ylab = "Index", xlab = "Times") # Code for generation Figure 14 acfret <- acf(diff(logdayprice)) acfret2 <- acf(diff(logdayprice)^2) par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3,2)) plot(acfret$acf[(2:11)], cex.main =.8, type = "h", main = "ACF Log ret", ylab = "ACF") plot(acfret2$acf[(2:11)], cex.main =.8, type = "h", main = "ACF squared Log ret", ylab = "ACF") ## Chunk 2: Arch Effect ## library("FinTS") # The package FinTS is used for detecting the Arch Effect logRetDay <- diff(logdayprice) ArchTest(logRetDay, lags = 5) ## Chunk 3: In and Out of Sample Data ## dayout <- 5 daysout <- obsinday * dayout Dataout <- allObs[(dimData - daysout + 1):dimData] Datain <- allObs[1:(dimData - daysout)] ## Chunk 4: Model Definition ## Cogparam11 <- list(a1 = 1.06e-02, b1 = 1.60e-02, a0 = 0.01, lambda = 1, alpha = sqrt(2), beta = 0, mu = 0, y01 = 0) Cogmodel11 <- setCogarch(p = 1, q = 1, measure = list(df = "rvgamma(z, lambda, alpha, beta, mu)"), measure.type = "code", XinExpr = TRUE) ## Chunk 5: In and Out of Sample Data comparison Fig. 15 ## incr.dataI <- diff(Datain) - mean(diff(Datain)) dataI <- setData(as.matrix(cumsum(c(0, incr.dataI))), delta = 1 / obsinday) incr.dataO <- diff(Dataout) - mean(diff(Dataout)) dataO <- setData(as.matrix(cumsum(c(0, incr.dataO))), delta = 1 / obsinday) # Code for generation Figure 15 par(mfrow = c(2, 1), mar = c(4.5, 4.5, 2.5, 2)) plot(dataI, main = "In-Sample Data") plot(dataO, main = "Out-Of-Sample Data") ## Chunk 6: Estimation COGARCH using Out-of-sample Data ## CogYuima11 <- setYuima(data = dataI, model = Cogmodel11) resCog11 <- gmm(yuima = CogYuima11, start = Cogparam11, Est.Incr = "IncrPar") summary(resCog11) TestStat <- Diagnostic.Cogarch(resCog11) ## Chunk 7: Estimated Increments Fig. 16 ## par(mfrow = c(1, 1), mar = c(4.5, 4.5, 3, 2)) plot(resCog11, main = "Increments") ## Chunk 8: Reconstruction of COGARCH trajectory with estimated increments Fig. 17 ## simCog <- simulate(resCog11) # Code for generation Figure 17 par(mfrow = c(1, 1), mar = c(4.5, 4.5, 3, 2)) plot(simCog, main = "Reconstructed COGARCH(1,1) sample path") ## Chunk 9: Arch effect on Innovations ## Incr.Lall <- resCog11@Incr.Lev Dates <- seq(1, length(Incr.Lall), by = obsinday) dayIncr <- diff(cumsum(Incr.Lall)[Dates]) ArchTest(dayIncr, lags = 5) ## Chunk 10: Simulated COGARCH trajectories using estimated Increments Fig. 18 ## mod <- resCog11@yuima@model Incr.L <- Incr.Lall[(length(Incr.Lall) - daysout):length(Incr.Lall)] param <- c(coef(resCog11), y01 = TestStat$meanStateVariable) set.seed(2) samp <- setSampling(Initial = 0, Terminal = dayout, n = obsinday * dayout) nrip <- 1000 AllGparam <- matrix(NA, obsinday * dayout, nrip) AllGboot <- matrix(NA, obsinday * dayout, nrip) for (i in 1:nrip) { pos <- as.integer(runif(obsinday * dayout, min = 1, max = length(Incr.L))) trajboot <- simulate(mod, true.parameter = param, Incr.L = as.matrix(Incr.L)[pos, ], samp = samp) traj <- simulate(mod,true.parameter = param, samp = samp) AllGparam[, i] <- traj@data@original.data[, 1] AllGboot[, i] <- trajboot@data@original.data[, 1] } # Code for generation Figure 18 par(mfrow = c(1, 2), mar = c(4.5, 4.5, 2.5, 2)) dtimeout <- index(dataO@original.data) / obsinday dataout <- dataO@original.data matplot(x = dtimeout, y = AllGparam, type = "l", ylim = c(-0.15, 0.05), ylab = "paths", xlab = "Days", main = "VG-COGARCH(1,1)", cex.main = 0.8) lines(x = dtimeout, y = dataout, type = "l",lwd = 3) matplot(x = dtimeout, y = AllGboot, type = "l", ylab = "paths",main = "NON-PAR COGARCH(1,1)", ylim = c(-0.15, 0.05), xlab = "Days", cex.main = 0.8) lines(x = dtimeout, y = dataout, type = "l", lwd = 3) ## Chunk 11: Estimation of a GARCH(1,1) on In sample data ## library("rugarch") daydatain <- diff(Datain[Dates]) spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0), include.mean = FALSE)) resGarch <- ugarchfit(data = daydatain, spec = spec) coef(resGarch) ## Chunk 12: Simulated GARCH trajectories using estimated Increments ## set.seed(2) bootp <- ugarchboot(resGarch, method = c("Partial", "Full")[1], n.ahead = 5, n.bootpred = 1000) series <- apply(bootp@fseries, 1, cumsum) ## Invisible Chunk 1: Code for Fig. 19 ## layout(matrix(c(1, 2, 3, 4, 5, 5), ncol = 2, byrow = TRUE),heights = rep.int(5, 2)) par(mai = rep(1, 4), mar = c(4.5, 4.5, 2.5, 2)) ncl<-30 min1 <- min(c(series[1, ], AllGboot[obsinday, ])) max1 <- max(c(series[1, ], AllGboot[obsinday, ])) break1 <- seq(min1, max1, length.out = ncl) GarchDens1 <- hist(series[1, ], plot = FALSE, breaks = break1) CogarchDens1 <- hist(AllGboot[obsinday, ], plot = FALSE, breaks = GarchDens1$breaks) plot(x = CogarchDens1$breaks[1:(length(CogarchDens1$count) + 1)], y = c(0, cumsum(CogarchDens1$count) / sum(CogarchDens1$count)), col = "blue", type = "l", xlab = "values", ylab = "CDF", xlim = c(-.021, .021), main="1 day CDF") lines(x = GarchDens1$breaks[1:(length(GarchDens1$count) + 1)], y = c(0, cumsum(GarchDens1$count) / sum(GarchDens1$count)), col = "red", type = "l") min5 <- min(c(series[5, ], AllGboot[obsinday * 5, ])) max5 <- max(c(series[5, ], AllGboot[obsinday * 5, ])) break5 <- seq(min5, max5, length.out = ncl) GarchDens5 <- hist(series[5, ], plot = FALSE, breaks = break5) CogarchDens5 <- hist(AllGboot[obsinday * 5,], plot = FALSE, breaks = break5) plot(x = CogarchDens5$breaks[1:(length(CogarchDens5$count) + 1)], y = c(0, cumsum(CogarchDens5$count) / sum(CogarchDens5$count)), col = "blue", type = "l", xlim = c(-.052, .052), xlab = "values", ylab = "CDF", main = "5 days CDF") lines(x = GarchDens5$breaks[1:(length(GarchDens5$count) + 1)], y = c(0, cumsum(GarchDens5$count) / sum(GarchDens5$count)), col = "red", type = "l") plot(x = CogarchDens1$mids, y = CogarchDens1$density, col = "blue", xlab = "values", type = "h", ylab = "PDF", xlim = c(-.021, .021), main = "1 day PDF", ylim = c(0, 85)) lines(x = GarchDens1$mids + 0.0003, y = GarchDens1$density, col = "red", type = "h") plot(x = CogarchDens5$mids,y = CogarchDens5$density, col = "blue", xlab = "values", type = "h", ylab = "PDF", xlim = c(-.052, .052), main = "5 days PDF", ylim = c(0, 30)) GarchDens5 <- hist(series[5, ], breaks = CogarchDens5$breaks, plot = FALSE) lines(x = GarchDens5$mids + 0.001, y = GarchDens5$density, col = "red", type = "h") par(mai = c(0, 0, 0, 0)) plot.new() legend("center", legend = c("COGARCH(1,1)", "GARCH(1,1)"), lty = 2, col = c("blue", "red"), bty = "n", cex = .75, horiz = TRUE) ## Invisible Chunk 2: Code for Fig. 20 ## par(mfrow = c(1, 2), mar = c(4.5, 4.5, 2.5, 2)) qqplot(x=series[1,],y=AllGboot[obsinday,], xlab="GARCH(1,1)", ylab="COGARCH(1,1)", main = "1 day QQ-PLOT", cex.main = 0.8) lines(x=break1,y=break1, col = "red") qqplot(x=series[5,],y=AllGboot[obsinday*5,], xlab="GARCH(1,1)", ylab="COGARCH(1,1)", main = "5 day QQ-PLOT", cex.main = 0.8) lines(x=break5,y=break5, col = "red") ## Invisible Chunk 3: Code for the Table before Section Conclusions ## mass <- NULL quan95 <- NULL quan90 <- NULL meanabs <- NULL quant10 <- NULL quant5 <- NULL for (i in 1:5) { min1 <- min(c(series[i, ], AllGboot[obsinday*i, ])) max1 <- max(c(series[i, ], AllGboot[obsinday*i, ])) break1 <- seq(min1, max1, length.out = ncl) GarchDens1 <- hist(series[i, ], plot = FALSE, breaks = break1) CogarchDens1 <- hist(AllGboot[obsinday*i, ], plot = FALSE, breaks = GarchDens1$breaks) GarchCDF1 <- cumsum(GarchDens1$count) / sum(GarchDens1$count) CogCDF1 <- cumsum(CogarchDens1$count) / sum(CogarchDens1$count) mass <- c(mass, max(abs(GarchCDF1 - CogCDF1))) quan95 <- c(quan95, quantile(abs(GarchCDF1 - CogCDF1), 0.95)) quan90 <- c(quan90, quantile(abs(GarchCDF1 - CogCDF1), 0.90)) meanabs <- c(meanabs, mean(abs(GarchCDF1 - CogCDF1))) quant10 <- c(quant10, quantile(abs(GarchCDF1 - CogCDF1), 0.10)) quant5 <- c(quant5, quantile(abs(GarchCDF1-CogCDF1),0.05)) } data <- matrix(0, 6, 5) data[1, ] <- mass data[2, ] <- quan95 data[3, ] <- quan90 data[4, ] <- meanabs data[5, ] <- quant10 data[6, ] <- quant5 colnames(data) <- c("I day", "II day", "III day", "IV day", "V day" ) rownames(data) <- c("max", "quantile.95", "quantile.90", "mean", "quantile.10", "quantile.05") tablequ <- as.table(data) show(tablequ) ################################################################# sessionInfo()