R/rgarch-armafor.R

In rgarch: Flexible GARCH modelling in R

#################################################################################
##
##   R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
##   This file is part of the R package rgarch.
##
##   The R package rgarch is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package rgarch 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 General Public License for more details.
##
#################################################################################
# forecast from an arma/garch-in-mean process
armaf = function(data, h, armaOrder, mu, ar, ma, inmean, im, epsx, include.mex, mxreg, mxf, n.ahead)
{
	# mxf is the extended external regressor matrix to
	# include their n.ahead forecasts provided in the
	# forecasting routine stage:
	# mxf: rbind(mexdata, external.forecasts$mregfor)
	N = length(data)
	x = c(data, rep(0,n.ahead))
	armap = armaOrder[1]
	armaq = armaOrder[2]
	if(include.mex){
		mu = mu + mxf%*%t(matrix(mxreg, ncol=dim(mxf)[2]))
	}
	# comment: we allow the arch-in-mean to decay (towards its unconditional value) 
	# at the rate of the variance forecast (->1 - persistence) rather than just using 
	# E[h]*inmean...
	# also note that h here is sigma not the variance
	mu = mu + inmean*(h^im)
	# correction for presence of xreg ( 20 Apr.2011 )
	for(i in 1:n.ahead){
		x[N+i] = mu[N+i] + ifelse(armap>0, sum(ar*(x[N+i-(1:armap)] - mu[N+i-(1:armap)])),0)
		for(j in 1:armaq){
			if(i-j>0){
				s = 0
			} else{
				s = ma[j]*epsx[N+i-j]
			}
			x[N+i] = x[N+i] + s
		}
	}
	return(x)
}

# forecast with arfima
arfimaf = function(data, h, armaOrder, mu, ar, ma, darfima, inmean, im, epsx, include.mex, mxreg, mxf, n.ahead)
{
	if(include.mex){
		mu = mu + mxf%*%t(matrix(mxreg, ncol = dim(mxf)[2]))
	}
	mu = mu + inmean*(h^im)
	x = c(data, rep(0, n.ahead))
	# demean data (includes regressors and arch-in-mean)
	xdata = x - mu
	
	N = length(data)
	armap = armaOrder[1]
	armaq = armaOrder[2]
	# generate fractional series
	zrf = .fracdiff(c( 1, rep(0, N + n.ahead - 1) ), darfima = darfima)
	if(armap>0){
		zar = rep(0, armap)
		for(i in 1:armap){
			zar[i] = -sum(rev(zrf[2:(N - i + 1)])*xdata[1:(N-i)])
		}
		# reverse back
		zar = rev(zar)
		# lagged data adjusted for fractional integration (i.e. subtracted
		# fractional mean)
		xzr = xdata[1:N] - c(rep(0, N-armap), zar)
	} else{
		# no need in the absence of armap since we have no lagged x terms
		# i.e. no 0 * x-mu = 0
		xzr = xdata[1:N]
	}
	for(i in 1:n.ahead){
		# add fractional mean back
		xzy = -sum(rev(zrf[2:(N + i)])*xdata[1:(N+i-1)])
		x[N+i] = (mu[N+i] + xzy) + ifelse(armap>0, sum(ar*(xzr[N+i-(1:armap)])), 0)
		for(j in 1:armaq){
			if(i-j>0){
				s = 0
			} else{
				s = ma[j]*epsx[N+i-j]
			}
			x[N+i] = x[N+i] + s
		}
		# forward adjust for the mean
		xdata[N+i] = x[N+i] - mu[N+i]
		# forward adjust for the lagged data (irrelevant when armap=0)
		xzr[N+i] = x[N+i] - mu[N+i] - xzy
	}
	return(x)
}