Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# METHOD: PREDICTION:
# predict.fGARCH Forecasts from an object of class 'fGARCH'
################################################################################
setMethod(f = "predict", signature(object = "fGARCH"), definition =
function(object, n.ahead = 10, trace = FALSE,
mse = c("cond","uncond"),
plot=FALSE, nx=NULL, crit_val=NULL, conf=NULL, ...,
p_loss = NULL # GNB: for ES and VaR
)
{
# A function implemented by Diethelm Wuertz
# Description:
# Prediction method for an object of class fGARCH
# Arguments:
# object an object of class fGARCH as returned by the
# function garchFit().
# n.ahead number of steps to be forecasted, an integer
# value, by default 10)
# trace should the prediction be traced? A logical value,
# by default FALSE)
# mse should the mean squared errors be conditional or unconditional
# plot should the predictions be plotted
# nx The number of observations to be plotted with the predictions
# (If plot is TRUE, the default value of nx is the sample
# size times 0.25.)
# crit_va If you want to set manually the critical values for
# the confidence intervals
# conf The confidence level for computing the critical values
# of the confidence intervals
# FUNCTION:
mse <- match.arg(mse)
# Retrieve "fit" from Parameter Estimation:
fit = object@fit
# Get ARMA(u,v)-GARCH(p,q) Order:
u = fit$series$order[1]
v = fit$series$order[2]
p = fit$series$order[3]
q = fit$series$order[4]
max.order = max(u, v, p, q)
# Get Start Conditions:
h.start = fit$series$h.start
llh.start = fit$series$llh.start
index = fit$params$index
params = fit$params$params
par = fit$par
Names = names(index)
for (Name in Names) params[Name] = par[Name]
Names = names(params)
# Retrieve From Initialized Parameters:
cond.dist = fit$params$cond.dist
# Extract the Parameters by Name:
leverage = fit$params$leverage
mu = params["mu"]
if (u > 0) {
ar = params[substr(Names, 1, 2) == "ar"]
} else {
ar = c(ar1 = 0)
}
if (v > 0) {
ma = params[substr(Names, 1, 2) == "ma"]
} else {
ma = c(ma1 = 0)
}
omega = params["omega"]
if (p > 0) {
alpha = params[substr(Names, 1, 5) == "alpha"]
} else {
alpha = c(alpha1 = 0)
}
if (p > 0 & leverage) {
gamma = params[substr(Names, 1, 5) == "gamma"]
} else {
gamma = c(gamma1 = 0)
}
if (q > 0) {
beta = params[substr(Names, 1, 4) == "beta"]
} else {
beta = c(beta1 = 0)
}
delta = params["delta"]
skew = params["skew"]
shape = params["shape"]
# Trace Parameters:
if (trace) {
cat("\nModel Parameters:\n")
print(c(mu, ar, ma, omega, alpha, gamma, beta, delta, skew, shape))
}
# Retrieve Series Lengths:
M = n.ahead
N = length(object@data)
# Get and Extend Series:
x = c(object@data, rep(mu, M))
h = c(object@h.t, rep(0, M))
z = c(fit$series$z, rep(mu, M))
# Forecast and Optionally Trace Variance Model:
var.model = fit$series$model[2]
# Forecast GARCH Variance:
if (var.model == "garch") {
if (trace) cat("\nForecast GARCH Variance:\n")
for (i in 1:M) {
h[N+i] = omega + sum(beta*h[N+i-(1:q)])
for (j in 1:p) {
if (i-j > 0) {
s = h[N + i - j]
} else {
s = z[N + i - j]^2
}
h[N+i] = h[N+i] + alpha[j] * s
}
}
}
# Forecast APARCH Variance:
if (var.model == "aparch") {
if (trace) cat("\nForecast APARCH Variance:\n")
for (i in 1:M) {
h[N+i] = omega + sum(beta*h[N+i-(1:q)])
for (j in 1:p) {
## 2024-01-30 GNB: TODO:
## it seems that kappa doesn't depend on i;
## so, kappa[1], ..., kappa[p] can be computed outside the i-loop
## and the formulas below use kappa[j]
kappa = garchKappa(cond.dist = cond.dist, gamma = gamma[j],
delta = delta, skew = skew, shape = shape)
if (i-j > 0) {
s = kappa * h[N + i - j]
} else {
s = (abs(z[N + i - j]) - gamma[j]*z[N + i - j])^delta
}
h[N+i] = h[N+i] + alpha[j] * s
}
}
}
# Forecast and Optionally Trace Mean Model:
# Note we set maxit=0 to get an object of class Arima with fixed
# init parameters ...
mu <- mu/(1-sum(ar))
ARMA <- arima(x = object@data, order = c(max(u, 1), 0, max(v, 1)),
init = c(ar, ma, mu), transform.pars = FALSE,
optim.control = list(maxit = 0))
prediction = predict(ARMA, n.ahead)
meanForecast = as.vector(prediction$pred)
if(mse=="uncond") {
meanError = as.vector(prediction$se)
} else {
# coefficients of h(t+1)
a_vec <- rep(0,(n.ahead))
hhat <- h[-(1:N)]^(2/delta[[1]]) #-> [[1]] to omit name of delta
u2 <- length(ar)
meanError <- hhat[1]
a_vec[1] = ar[1] + ma[1]
meanError <- na.omit(c(meanError,sum(hhat[1:2]*c(a_vec[1]^2,1))))
if ((n.ahead - 1) > 1) {
for( i in 2:(n.ahead - 1)) {
a_vec[i] <- ar[1:min(u2,i-1)]*a_vec[(i-1):(i-u2)] +
ifelse(i>u,0,ar[i]) + ifelse(i>v,0,ma[i])
meanError <- na.omit(c(meanError,
sum(hhat[1:(i+1)]*c(a_vec[i:1]^2,1))))
}
}
meanError <- sqrt(meanError)
}
if (trace) {
cat("\nForecast ARMA Mean:\n")
print(ARMA)
cat("\n")
print(prediction)
}
# Standard Deviations:
standardDeviation = h^(1/delta)
# Plotting the predictions
if (plot) {
if(is.null(nx))
nx <- round(length(object@data)*.25)
t <- length(object@data)
x <- c(object@data[(t-nx+1):t],meanForecast)
# Computing the appropriate critical values
if (is.null(conf))
conf <- 0.95
if (is.null(crit_val)) {
if (object@fit$params$cond.dist=="norm") {
crit_valu <- qnorm(1-(1-conf)/2)
crit_vald <- qnorm((1-conf)/2)
}
if (object@fit$params$cond.dist=="snorm") {
crit_valu <- qsnorm(1-(1-conf)/2,xi=coef(object)["skew"])
crit_vald <- qsnorm((1-conf)/2,xi=coef(object)["skew"])
}
if (object@fit$params$cond.dist=="ged") {
crit_valu <- qged(1-(1-conf)/2,nu=coef(object)["shape"])
crit_vald <- qged((1-conf)/2,nu=coef(object)["shape"])
}
if (object@fit$params$cond.dist=="sged") {
crit_valu <- qsged(1-(1-conf)/2,nu=coef(object)["shape"],
xi=coef(object)["skew"])
crit_vald <- qsged((1-conf)/2,nu=coef(object)["shape"],
xi=coef(object)["skew"])
}
if (object@fit$params$cond.dist=="std") {
crit_valu <- qstd(1-(1-conf)/2,nu=coef(object)["shape"])
crit_vald <- qstd((1-conf)/2,nu=coef(object)["shape"])
}
if (object@fit$params$cond.dist=="sstd") {
crit_valu <- qsstd(1-(1-conf)/2,nu=coef(object)["shape"],
xi=coef(object)["skew"])
crit_vald <- qsstd((1-conf)/2,nu=coef(object)["shape"],
xi=coef(object)["skew"])
}
if (object@fit$params$cond.dist=="snig") {
crit_valu <- qsnig(1-(1-conf)/2,zeta=coef(object)["shape"],
rho=coef(object)["skew"])
crit_vald <- qsnig((1-conf)/2,zeta=coef(object)["shape"],
rho=coef(object)["skew"])
}
if (object@fit$params$cond.dist=="QMLE") {
e <- sort(object@residuals/object@sigma.t)
crit_valu <- e[round(t*(1-(1-conf)/2))]
crit_vald <- e[round(t*(1-conf)/2)]
}
} else {
if (length(crit_val)==2) {
crit_valu <- crit_val[2]
crit_vald <- crit_val[1]
}
if (length(crit_val)==1) {
crit_valu <- abs(crit_val)
crit_vald <- -abs(crit_val)
}
}
int_l <- meanForecast+crit_vald*meanError
int_u <- meanForecast+crit_valu*meanError
ylim_l <- min(c(x,int_l)*(.95))
ylim_u <- max(c(x,int_u)*(1.05))
plot(x,type='l',ylim=c(ylim_l,ylim_u))
title("Prediction with confidence intervals")
lines((nx+1):(nx+n.ahead), meanForecast, col = 2, lwd = 2)
lines((nx+1):(nx+n.ahead), int_l, col = 3, lwd = 2)
lines((nx+1):(nx+n.ahead), int_u, col = 4, lwd = 2)
polygon(c((nx+1):(nx+n.ahead),(nx+n.ahead):(nx+1)),
c(int_l, int_u[n.ahead:1]),
border = NA, density = 20, col = 5, angle = 90)
es1 <- as.expression(substitute(hat(X)[t+h] + crit_valu*sqrt(MSE),
list(crit_valu=round(crit_valu,3))))
es2 <- as.expression(substitute(hat(X)[t+h] - crit_vald*sqrt(MSE),
list(crit_vald=abs(round(crit_vald,3)))))
es3 <- expression(hat(X)[t+h])
legend("bottomleft",c(es3,es2,es1),col=2:4,lty=rep(1,3),lwd=rep(2,3))
grid()
}
## Result:
forecast <- data.frame(
meanForecast = meanForecast,
meanError = meanError,
standardDeviation = standardDeviation[-(1:N)])
if(plot)
forecast = data.frame(
forecast,
lowerInterval = int_l,
upperInterval = int_u)
## if(plot) {
## forecast = data.frame(
## meanForecast = meanForecast,
## meanError = meanError,
## standardDeviation = standardDeviation[-(1:N)],
## lowerInterval = int_l,
## upperInterval = int_u)
## } else {
## forecast = data.frame(
## meanForecast = meanForecast,
## meanError = meanError,
## standardDeviation = standardDeviation[-(1:N)])
## }
## 2024-01-31 GNB: VaR and ES - experimental
if(!is.null(p_loss)) {
cond_dist <- object@fit$params$cond.dist
mu_t <- meanForecast
sigma_t <- standardDeviation[-(1:N)]
qf <- qfun_fGarch(cond_dist, coef(object))
predVaR <- cvar::VaR_qf(qf, p_loss, intercept = mu_t, slope = sigma_t)
predES <- cvar::ES(qf, p_loss, intercept = mu_t, slope = sigma_t)
forecast <- data.frame(forecast, VaR = predVaR, ES = predES)
attr(forecast, "p_loss") <- p_loss
}
# Return Value:
forecast
})
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.