################################################### ### Load the package and produce two graphical ### diagnostic plots. ################################################### library("evmix") data("spto87", package = "evir") mrlplot(spto87, try.thresh = c(2.6, 1.5, 0.4)) tshapeplot(spto87, try.thresh = c(2.7, 1.8, 0.4), ylim = c(-1, 1), legend = "bottomleft") pickandsplot(spto87, try.thresh = c(1.8, 1.3, 1.05), ylim = c(-1, 1)) ################################################### ### Plot example of normal+GPD probability density ### function and cumulative distribution function. ################################################### x <- seq(-5, 5, 0.001) y <- dnormgpd(x, u = 1.5, sigmau = 0.4) plot(x, y, xlab = "x", ylab = "Density f(x)", type = "l") abline(v = 1.5, lty = 2) y <- pnormgpd(x, u = 1.5, sigmau = 0.4) plot(x, y, xlab = "x", ylab = "Distribution F(x)", type = "l") abline(v = 1.5, lty = 2) ################################################### ### Simulate 50 data points from standard normal, ### estimate bandwidth from standard KDE and ### overlay KDE over a density histogram. ################################################### set.seed(0) xdata <- rnorm(50) fit <- fkden(xdata) xp <- seq(-4, 4, 0.01) hist(xdata, 10, freq = FALSE, xlim = c(-4, 4), ylim = c(0, 0.7)) rug(xdata) lines(xp, dnorm(xp), col = "black") with(fit, lines(xp, dkden(xp, x, lambda), col = "red")) lines(density(xdata), lty = 2, col = "green") legend("topright", c("True Density", "KDE fitted by MLE in evmix", "KDE using density, default bw"), lty = c(1, 1, 2), col = c("black", "red", "green")) ################################################### ### Simulate 100 datapoints from Student-t(df = 4), ### estimate bandwidth for standard KDE from entire ### dataset, then use a trimmed estimator of the ### bandwidth. ### Overlay the two fits over the density histogram. ################################################### set.seed(0) xdata <- rt(100, df = 2) fit <- fkden(xdata) intails <- ((xdata <= -3) | (xdata >= 3)) fit.notails <- fkden(xdata[!intails], extracentres = xdata[intails]) xp <- seq(-10, 10, 0.01) hist(xdata, seq(floor(min(xdata)), ceiling(max(xdata)), 0.5), freq = FALSE, xlim = c(-10, 10)) rug(xdata) lines(xp, dt(xp , df = 2), col = "black") with(fit, lines(xp, dkden(xp, x, lambda), col = "blue")) with(fit.notails, lines(xp, dkden(xp, x, lambda), col = "red")) legend("topright", c("True Density", "KDE fitted by MLE", "KDE if tails trimmed"), col = c("black", "blue", "red")) ################################################### ### Simulate 500 data points from a gamma(1,2), ### estimate the bandwidth using the "reflection" ### boundary corrected KDE method. ### Use this estimate bandwidth to estimate the ### the density using the reflection, ### renomalization and simple boundary correction ### methods. ################################################### set.seed(0) xdata <- rgamma(500, shape = 1, scale = 2) fit <- fbckden(xdata, linit = 1, bcmethod = "reflect") xp <- seq(-1, 10, 0.01) hist(xdata, 20, freq = FALSE, ylim = c(0, 0.6)) rug(xdata) lines(xp, dgamma(xp, shape = 1, scale = 2), col = "black") with(fit, lines(xp, dbckden(xp, xdata, lambda, bcmethod = bcmethod), col = "red")) with(fit, lines(xp, dbckden(xp, xdata, lambda, bcmethod = "renorm"), col = "brown")) with(fit, lines(xp, dbckden(xp, xdata, lambda, bcmethod = "simple"), col = "blue")) lines(density(xdata), lty = 2, col = "green") legend("topright", c("True Density", "BC KDE using reflect", "BC KDE using renorm", "BC KDE using simple", "KDE Using density"), lty = c(rep(1, 4), 2), col = c("black", "red", "brown", "blue", "green")) ################################################### ### Load the Standard and Poor's daily returns data. ### Create a time series plot of the data. ################################################### data("spto87", package = "evir") plot(attr(spto87, "times"), spto87, type = "l", xlab = "Date", ylab = "Daily log Returns") ################################################### ### Fit the normal+GPD model, with bulk model and ### parameterized tail fraction approach. ################################################### fit.bulk <- fnormgpd(spto87, useq = seq(0.2, 2, 0.1), fixedu = TRUE) fit.para <- fnormgpd(spto87, useq = seq(0.2, 2, 0.1), fixedu = TRUE, phiu = FALSE) ################################################### ### Overlay the two fitted normal+GPD models over ### the density histogram. ################################################### hist(spto87, freq = FALSE, breaks = 50, xlab = "Daily log Returns", xlim = c(-4, 4), ylim = c(0, 0.8)) xp <- seq(-7, 7, 0.01) lines(xp, dnorm(xp, mean(spto87), sd(spto87)), col = "red") fx.bulk <- with(fit.bulk, dnormgpd(xp, nmean, nsd, u, sigmau, xi)) lines(xp, fx.bulk, col = "blue", lwd = 2) abline(v = fit.bulk$u, col = "blue", lwd = 2) fx.para <- with(fit.para, dnormgpd(xp, nmean, nsd, u, sigmau, xi, phiu)) lines(xp, fx.para, col = "green", lty = 2, lwd = 2) abline(v = fit.para$u, col = "green", lty = 2, lwd = 2) ################################################### ### Produce the usual return level plot, to explore ### the model fit in the upper tail. ################################################### rlplot(fit.bulk) rlplot(fit.para) ################################################### ### Focus on a key region of the return level plot. ################################################### rlplot(fit.bulk, rplim = c(10, 300), rllim = c(1, 3)) rlplot(fit.para, rplim = c(10, 300), rllim = c(1, 3)) ################################################### ### Show the sensitivity of the optimization ### algorithm to the initial value of the threshold. ################################################### fit.default <- fnormgpd(spto87) c(fit.default$u, quantile(spto87, 0.9)) ################################################### ### Fit a two tailed GPD+normal+GPD model to the ### Standard and Poor daily returns. ### Overlay the fit on the density histogram. ################################################### fit.gng <- fgng(spto87, ulseq = seq(-2, -0.2, 0.1), urseq = seq(0.2, 2, 0.1), fixedu = TRUE) 2 * (fit.bulk$nllh - fit.gng$nllh) qchisq(0.95, df = 3) hist(spto87, freq = FALSE, breaks = 50, xlab = "Daily Closing Returns", xlim = c(-4, 4), ylim = c(0, 0.8)) fx <- with(fit.gng, dgng(xp, nmean, nsd, ul, sigmaul, xil, phiul = TRUE, ur, sigmaur, xir, phiur = TRUE)) lines(xp, fx) abline(v = c(fit.gng$ul, fit.gng$ur)) ################################################### ### Produce the usual model diagnostic quartet. ################################################### evmix.diag(fit.gng) ################################################### ### Fit a two tailed GPD+KDE+GPD model to the ### Standard and Poor daily returns. ### Overlay the fit on the density histogram. ################################################### fit.gkg <- fgkg(spto87, ulseq = fit.gng$ul, urseq = fit.gng$ur, fixedu = TRUE) hist(spto87, freq = FALSE, breaks = 50, xlab = "Daily Closing Returns", xlim = c(-4, 4), ylim = c(0, 0.8)) fx <- with(fit.gkg, dgkg(xp, spto87, lambda, ul, sigmaul, xil, phiul = TRUE, ur, sigmaur, xir, phiur = TRUE)) lines(xp, fx) abline(v = c(fit.gkg$ul, fit.gkg$ur)) ################################################### ### Explore the upper tail fit using the return ### level plot. ################################################### rlplot(fit.gkg, alpha = NULL, rplim = c(1, 11000), rllim = c(0, 5)) ################################################### ### Reproduce Table 1 ################################################### library("evmix") data("spto87", package = "evir") # Fit GPD to range of thresholds gpd04 = fgpd(spto87, u = 0.4) gpd10 = fgpd(spto87, u = 1) gpd13 = fgpd(spto87, u = 1.3) gpd15 = fgpd(spto87, u = 1.5) gpd18 = fgpd(spto87, u = 1.8) gpd26 = fgpd(spto87, u = 2.6) gpd27 = fgpd(spto87, u = 2.7) # Compile table 1 table1results = data.frame(u1 = gpd04$u, xi1 = gpd04$xi, xi1se = gpd04$se[2], u2 = gpd15$u, xi2 = gpd15$xi, xi2se = gpd15$se[2], u3 = gpd26$u, xi3 = gpd26$xi, xi3se = gpd26$se[2]) table1results = rbind(table1results, c(with(gpd04, c(u, xi, se[2])), with(gpd18, c(u, xi, se[2])), with(gpd27, c(u, xi, se[2])))) table1results = rbind(table1results, c(with(gpd10, c(u, xi, se[2])), with(gpd13, c(u, xi, se[2])), with(gpd18, c(u, xi, se[2])))) print(table1results) ################################################### ### Reproduce Table 2 ################################################### # Load Standard and Poors 1987 data from evir package library("evmix") data("spto87", package = "evir") # Fit normal+GPD mixture bulktime = system.time(fit.bulk <- fnormgpd(spto87, useq = seq(0.2, 2, 0.1), fixedu = TRUE)) paratime = system.time(fit.para <- fnormgpd(spto87, useq = seq(0.2, 2, 0.1), fixedu = TRUE, phiu = FALSE)) gngtime = system.time(fit.gng <- fgng(spto87, ulseq = seq(-2, -0.2, 0.1), urseq = seq(0.2, 2, 0.1), fixedu = TRUE)) # Compile table 2 table2results = with(fit.bulk, data.frame(ur = u, xir = xi, xirse = se[4], nllh = -nllh, npar = length(mle) + 1, computetime = bulktime[3])) # + 1 parameters for u table2results = rbind(table2results, c(with(fit.para, c(u, xi, se[4], -nllh, length(mle) + 2)), paratime[3])) # + 2 parameters for u and phiu table2results = rbind(table2results, c(with(fit.gng, c(ur, xir, se[6], -nllh, length(mle) + 2)), gngtime[3])) # + 2 parameters for ul and ur print(table2results)