Nothing
#################################################################################
##
## 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.
##
#################################################################################
# A set of functions to compute VARX models for internal use
# fit/filter, forecast and simulate.
# Partical Credit goes to Pfaff for his vars package, code from which has been
# partially copied (particularly first part of varxfilter).
# Because of the reduced functionality required of
# lm is not called. Also, varxsim is completely new, and the forecast routine implements
# rolling forecasts with out.sample data.
varxfilter = function(X, p, exogen = NULL, Bcoef = NULL, robust = FALSE, gamma = 0.25, delta = 0.01, nc = 10, ns = 500)
{
X = as.matrix(X)
if (any(is.na(X))) stop("\nvarxfit:-->error: NAs in X.\n")
if (ncol(X) < 2) stop("\nvarxfit:-->error: The matrix 'X' should contain at least two variables.\n")
if (is.null(colnames(X))) colnames(X) = paste("X", 1:ncol(X), sep = "")
colnames(X) = make.names(colnames(X))
obs = dim(X)[1]
K = dim(X)[2]
xsample = obs - p
Xlags = embed(X, dimension = p + 1)[, -(1:K)]
temp1 = NULL
for (i in 1:p) {
temp = paste(colnames(X), ".l", i, sep = "")
temp1 = c(temp1, temp)
}
colnames(Xlags) = temp1
Xend = X[-c(1:p), ]
rhs = cbind( Xlags, rep(1, xsample))
colnames(rhs) <- c(colnames(Xlags), "const")
if( !(is.null(exogen)) ) {
exogen = as.matrix(exogen)
if (!identical(nrow(exogen), nrow(X))) {
stop("\nvarxfit:-->error: Different row size of X and exogen.\n")
}
XK = dim(exogen)[2]
if (is.null(colnames(exogen))) colnames(exogen) = paste("exo", 1:ncol(exogen), sep = "")
colnames(exogen) = make.names(colnames(exogen))
tmp = colnames(rhs)
rhs = cbind(rhs, exogen[-c(1:p), ])
colnames(rhs) = c(tmp, colnames(exogen))
} else{
XK = 0
}
datamat = as.matrix(rhs)
colnames(datamat) = colnames(rhs)
if( is.null( Bcoef ) ){
if( robust ){
sol = robustvar(data = X, exogen = exogen, lags = p, alpha = gamma, ns = ns, nc = nc, delta = delta)
Bcoef = t(sol$betaR)
Bcoef = cbind(Bcoef[,2:(p*K+1)], Bcoef[,1], if(!is.null(exogen)) Bcoef[,(2+p*K):(1+p*K+XK)] else NULL)
colnames(Bcoef) = colnames(rhs)
rownames(Bcoef) = colnames(X)
xfitted = t( Bcoef %*% t( datamat ) )
xresiduals = tail(X, obs - p) - xfitted
sigma2 = diag( sol$sigmaR )
# I think we need to correct this to get the corrected values from the mlts function
xp = solve( t(datamat)%*%datamat, diag(p*K + XK + 1) )
Bcov = lapply( as.list(sigma2), FUN = function(x) x * xp )
# calculate standard errors and t-stats
se = sapply(Bcov, FUN = function(x) sqrt(diag(x)))
rownames(se) = colnames(rhs)
tstat = t(Bcoef) / se
pstat = 2*(1-pnorm(abs(tstat)))
} else{
Bcoef = matrix(NA, ncol = dim(datamat)[2], nrow = K)
Bcoef = t( apply(Xend, 2, FUN = function(x) solve(t(datamat) %*% datamat) %*% t(datamat) %*% x ) )
colnames(Bcoef) = colnames(rhs)
rownames(Bcoef) = colnames(X)
xfitted = t( Bcoef %*% t( datamat ) )
xresiduals = tail(X, obs - p) - xfitted
sigma2 = apply( xresiduals, 2, FUN = function(x) (t(x) %*% x)/(obs - p*K - XK - 1) )
xp = solve( t(datamat)%*%datamat, diag(p*K + XK + 1) )
Bcov = lapply( as.list(sigma2), FUN = function(x) x * xp )
# calculate standard errors and t-stats
se = sapply(Bcov, FUN = function(x) sqrt(diag(x)))
rownames(se) = colnames(rhs)
tstat = t(Bcoef) / se
pstat = 2*(1-pnorm(abs(tstat)))
}
} else{
colnames(Bcoef) = colnames(rhs)
rownames(Bcoef) = colnames(X)
xfitted = t( Bcoef %*% t( datamat ) )
xresiduals = tail(X, obs - p) - xfitted
sigma2 = apply( xresiduals, 2, FUN = function(x) (t(x) %*% x)/(obs - p*K - XK - 1) )
xp = solve( t(datamat)%*%datamat, diag(p*K + XK + 1) )
Bcov = lapply( as.list(sigma2), FUN = function(x) x * xp )
# calculate standard errors and t-stats
se = sapply(Bcov, FUN = function(x) sqrt(diag(x)))
rownames(se) = colnames(rhs)
tstat = t(Bcoef) / se
pstat = 2*(1-pnorm(abs(tstat)))
}
ans = list( Bcoef = Bcoef, xfitted = xfitted, xresiduals = xresiduals, Bcov = Bcov, se = se, tstat = tstat, pstat = pstat, lag = p,
mxn = XK )
return( ans )
}
# work in progress
.uncmeanvar = function(Bcoef, p, mxn)
{
# we work with constant:
m = dim(Bcoef)[1]
Id = diag(m)
C = Bcoef[, p * m + 1]
if(mxn>0){
EX = Bcoef[, ( p * m + 2 ):( p * m + mxn + 1)]
Z = Bcoef[,-( ( p*m + 1 ):( p*m + mxn + 1) )]
#for(i in 1:mxn) C = C + (EX[, i])
} else{
EX = NULL
Z = Bcoef[,-((p*m)+1)]
}
idx = t( .embed(1:(p*m), m, m, ascending = TRUE) )
for(i in 1:p){
Id = Id - Z[,idx[,i]]
}
(1/Id) %*% C
}
varxforecast = function(X, Bcoef, p, out.sample, n.ahead, n.roll, mregfor)
{
X = as.matrix(X)
X.orig = X
n.roll = n.roll + 1
m = ncol(X)
n = nrow(X)
X = X.orig[1:(n - out.sample), ,drop = F]
n = nrow(X)
if(!is.null(mregfor)){
mregfor = as.matrix(mregfor)
mxn = ncol(mregfor)
} else{
mxn = 0
}
if(n.roll == 1){
dmat1 = NULL
for(i in 0:(p-1))
{
dmat1 = cbind(dmat1, .lagx(X, n.lag = i, pad = 0) )
}
# for n.ahead > 1 and n.roll == 1
meanforc = matrix(NA, ncol = m, nrow = n.ahead)
dmat1 = as.numeric( tail(dmat1, 1) )
for(i in 1:n.ahead){
Z = cbind(matrix(dmat1[1 : (m * p)], nrow = 1), matrix(1, ncol = 1, nrow = 1), if(!is.null(mregfor)) matrix(as.numeric(mregfor[i, ]), nrow = 1 ) else NULL)
meanforc[i, ] = Bcoef %*% t(Z)
dmat1 = c(meanforc[i, ], dmat1)
}
} else{
dmat2 = NULL
N = dim(X.orig)[1] - out.sample
for(i in 0:(p-1))
{
dmat2 = cbind(dmat2, .lagx(X.orig, n.lag = i, pad = 0) )
}
# n.roll
dmat = dmat2[N, ]
meanforc = matrix(NA, ncol = m, nrow = n.roll)
for(i in 1:n.roll){
Z = cbind(matrix(dmat[1 : (m * p)], nrow = 1), matrix(1, ncol = 1, nrow = 1), if(!is.null(mregfor)) matrix(as.numeric(mregfor[i, ]), nrow = 1 ) else NULL)
meanforc[i, ] = Bcoef %*% t(Z)
if( i < n.roll) dmat = c(dmat2[N+i, ], dmat)
}
}
return( meanforc )
}
varxsim = function(X, Bcoef, p, n.sim, n.start, prereturns, resids, mexsimdata)
{
X = as.matrix(X)
n = dim(X)[1]
m = dim(X)[2]
ns = n.start + n.sim
if(!is.null(mexsimdata)){
mexsimdata = as.matrix(mexsimdata)
mxn = ncol(mexsimdata)
} else{
mxn = 0
Bcoef = Bcoef[,1:(m*p+1)]
}
dmat2 = NULL
if( is.null(prereturns) ){
for(i in 0:(p-1))
{
dmat2 = cbind(dmat2, .lagx(X, n.lag = i, pad = 0) )
}
} else{
for(i in 0:(p-1))
{
dmat2 = cbind(dmat2, .lagx(prereturns, n.lag = i, pad = 0) )
}
}
if( is.null(resids) ){
resids = matrix(0, ncol = m, nrow = ns)
}
dmat = tail(dmat2, 1)
meansim = matrix(NA, ncol = m, nrow = ns)
for(i in 1:ns){
Z = cbind(matrix(dmat[1 : (m * p)], nrow = 1), matrix(1, ncol = 1, nrow = 1), if(!is.null(mexsimdata)) matrix(as.numeric(mexsimdata[i, ]), nrow = 1 ) else NULL)
meansim[i, ] = Bcoef %*% t(Z) + resids[i, ]
if( i < ns ) dmat = c(meansim[i, ], dmat)
}
meansim = tail(meansim, n.sim)
rownames(meansim) = 1:n.sim
return( meansim )
}
varxsimXX = function(X, Bcoef, p, m.sim, prereturns, resids, mexsimdata)
{
X = as.matrix(X)
n = dim(X)[1]
m = dim(X)[2]
if(!is.null(mexsimdata)){
mexsimdata = as.matrix(mexsimdata)
mxn = ncol(mexsimdata)
# quick check to see if mxn = 1 (in which case it is inverted)
if( mxn == m.sim){
mexsimdata = t(mexsimdata)
mxn = ncol(mexsimdata)
}
mxBcoef = matrix(Bcoef[,(m*p+2):(m*p+1+mxn)], ncol = mxn)
Bcoef = Bcoef[,1:(m*p+1)]
} else{
mxn = 0
Bcoef = Bcoef[,1:(m*p+1)]
}
dmat2 = NULL
if( is.null(prereturns) ){
for(i in 0:(p-1))
{
dmat2 = cbind(dmat2, .lagx(X, n.lag = i, pad = 0) )
}
} else{
for(i in 0:(p-1))
{
dmat2 = cbind(dmat2, .lagx(prereturns, n.lag = i, pad = 0) )
}
}
if( is.null(resids) ){
resids = matrix(0, ncol = m, nrow = m.sim)
}
dmat = tail(dmat2, 1)
#meansim = matrix(NA, ncol = m, nrow = m.sim)
Z = cbind(matrix(dmat[1 : (m * p)], nrow = 1), matrix(1, ncol = 1, nrow = 1) )
meansim = matrix(Bcoef %*% t(Z), ncol = m, nrow = m.sim, byrow = TRUE)
if(!is.null(mexsimdata)){
for(i in 1:m.sim){
meansim[i,] = meansim[i,] + matrix(mxBcoef %*% mexsimdata[i, ], ncol = m, byrow = TRUE)
}
}
meansim = meansim + resids
return( meansim )
}
## VAR lag selection copied and amended from the 'vars' package
.varxselect = function (y, lag.max = 10, exogen = NULL)
{
y = as.matrix(y)
if (any(is.na(y))) stop("\nNAs in Data.\n")
colnames(y) = make.names(colnames(y))
K = ncol(y)
lag.max = abs(as.integer(lag.max))
lag = abs(as.integer(lag.max + 1))
ylagged = embed(y, lag)[, -c(1:K)]
yendog = y[-c(1:lag.max), ]
sample = nrow(ylagged)
rhs = rep(1, sample)
if (!(is.null(exogen))) {
exogen = as.matrix(exogen)
if (!identical(nrow(exogen), nrow(y))) { stop("\nDifferent row size of Data and exogenous regressors.\n") }
if (is.null(colnames(exogen))) {
colnames(exogen) = paste("exo", 1:ncol(exogen),
sep = "")
}
colnames(exogen) = make.names(colnames(exogen))
rhs = cbind(rhs, exogen[-c(1:lag.max), ])
}
idx = seq(K, K * lag.max, K)
detint = ncol(as.matrix(rhs))
criteria = matrix(NA, nrow = 4, ncol = lag.max)
rownames(criteria) = c("AIC", "HQ", "SC", "FPE")
colnames(criteria) = paste(seq(1:lag.max))
for (i in 1:lag.max) {
ys.lagged = cbind(ylagged[, c(1:idx[i])], rhs)
sampletot = nrow(y)
nstar = ncol(ys.lagged)
resids = resid(lm(yendog ~ -1 + ys.lagged))
sigma.det = det(crossprod(resids)/sample)
criteria[1, i] = log(sigma.det) + (2/sample) * (i *
K^2 + K * detint)
criteria[2, i] = log(sigma.det) + (2 * log(log(sample))/sample) *
(i * K^2 + K * detint)
criteria[3, i] = log(sigma.det) + (log(sample)/sample) *
(i * K^2 + K * detint)
criteria[4, i] = ((sample + nstar)/(sample - nstar))^K *
sigma.det
}
order = apply(criteria, 1, which.min)
return(list(selection = order, criteria = criteria))
}
.varcovres = function(X, mxn, p, resids){
n = dim(X)[1]
m = dim(X)[2]
# series x lags + constant + exogenous
nk = m * p + 1 + mxn
cov(resids) * (n - 1) / (n - nk)
}
################################################################################
# Multivariate Regression
# Input:
# x: data-matrix (n,p)
# y: data-matrix (n,q)
# gamma: proportion of trimming
# arguments:
# - ns : contains number of subsets; default=5000
# - nc : number of C-steps; default=10
# - delta : critical value for Reweighted estimator, deafult=0.01
# Output:
# beta : matrix (p,q) of MLTS-regression coefficients
# sigma: matrix (q,q) containing MLTS-residual covariance
# dres : residual distances (n,1) w.r.t. initial fit
# betaR : matrix (p,q) of RMLTS-regression coefficients
# sigmaR: matrix (q,q) containing RMLTS-residual covariance
# dresR : residual distances (n,1) w.r.t. RMLTS
# Remark: if intercept needed, add a column of ones to the x-matrix
#
# Ref: Agullo,J., Croux, C., and Van Aelst, S. (2008)
# The Multivariate Least Trimmed Squares Estimator,
# Journal of multivariate analysis
#
# Author: Kristel Joossens
#
#
###########
# EXAMPLE #
###########
# n=1000;
# p=5;
# q=2;
# x = cbind(rep(1,n),array(rnorm(n*(p-1)),dim=c(n,p-1)))
# beta = cbind(c(2,1,4,2,1),c(5,2,1,.5,2))
# y = x %*% beta +array(rnorm(n*q),dim=c(n,q))
# gamma = 0.25
# out = mlts(x,y,gamma)
################################################################################
mlts <- function(x,y,gamma,ns=500,nc=10,delta=0.01)
{
d <- dim(x); n <- d[1]; p <- d[2]
q <- ncol(y)
h <- floor(n*(1-gamma))+1
obj0 <- 1e10
for (i in 1:ns)
{ sorted <- sort(runif(n),na.last = NA,index.return=TRUE)
istart <- sorted$ix[1:(p+q)]
xstart <- x[istart,]
ystart <- y[istart,]
bstart <- solve(t(xstart)%*%xstart,t(xstart)%*%ystart)
sigmastart <- (t(ystart-xstart%*%bstart))%*%(ystart-xstart%*%bstart)/q
for (j in 1:nc)
{ res <- y - x %*% bstart
tres <- t(res)
dist2 <- colMeans(solve(sigmastart,tres)*tres)
sdist2 <- sort(dist2,na.last = NA,index.return = TRUE)
idist2 <- sdist2$ix[1:h]
xstart <- x[idist2,]
ystart <- y[idist2,]
bstart <- solve(t(xstart)%*%xstart,t(xstart)%*%ystart)
sigmastart <- (t(ystart-xstart%*%bstart))%*%(ystart-xstart%*%bstart)/(h-p)
}
obj <- det(sigmastart)
if (obj < obj0)
{ result.beta <- bstart
result.sigma <- sigmastart
obj0 <- obj
}
}
cgamma <- (1-gamma)/pchisq(qchisq(1-gamma,q),q+2)
result.sigma <- cgamma * result.sigma
res <- y - x %*% result.beta
tres<-t(res)
result.dres <- colSums(solve(result.sigma,tres)*tres)
result.dres <- sqrt(result.dres)
qdelta <- sqrt(qchisq(1-delta,q))
good <- (result.dres <= qdelta)
xgood <- x[good,]
ygood <- y[good,]
result.betaR <- solve(t(xgood)%*%xgood,t(xgood)%*%ygood)
result.sigmaR <- (t(ygood-xgood%*%result.betaR)) %*%
(ygood-xgood%*%result.betaR)/(sum(good)-p)
cdelta <- (1-delta)/pchisq(qdelta^2,q+2)
result.sigmaR<-cdelta*result.sigmaR
resR<-y-x%*%result.betaR
tresR<-t(resR)
result.dresR <- colSums(solve(result.sigmaR,tresR)*tresR)
result.dresR <- sqrt(result.dresR)
list(beta=result.beta,sigma=result.sigma,dres=result.dres,
betaR=result.betaR,sigmaR=result.sigmaR,dresR=result.dresR)
}
################################################################################
# Robust estimation of a VAR model using the robust MLST estimator.
# Robust lag length criteria for the give lag are reported.
# The user can as for the robust impulse response functions and their
# analytic confidence bounds.
#
# Input arguments:
# data : multivariate time series (T x nvar) (nvar=number of variables)
# lags : order of the VAR model
# The three possible extra arguments can be
# alpha : trimming portion to obtain the MLTS estimator
# the default value is 0.25
# delta : trimming portion to obtain the RMLTS estimator (based on MLTS)
# the default value is 0.01. If you want to redefine delta, you
# have to redefine alpha as previous argument
# nstep : number of steps: If are only interested in the effect of a unit
# shock at an innovation to observations not far/far in the
# future. You have to take nstep small(nstep=20)/large(nsetp=200)
# The combinations for the last 3 input arguments might be as follows
# nothing - alpha - alpha,delta - alpha,delta,nstep -
# nstep - nstep,alpha - nstep,alpha,delta
# eg. robustvar(data,lags,nstep,alpha)
# Output arguments:
# beta : the estimated coefficients
# AIC : The Akaike criterion
# HQ : The Hannan Quinn criterion
# SC : The Schwarz criterion
# If nstep is not given, no IRF will be created.
# If nstep is given: figure 1 represent the IRFs and their Analytic CB for
# VAR(lags)-model using OLS.
# Z : the estimated IRF
# UB : The 95# upper bound for the IRF
# LB : The 95# lower bound for the IRF
# robdis : contains the robust Mahalanobis distances
robustvar = function(data, exogen = NULL, lags = 2, alpha = 0.01, ns = 500, nc = 10, delta = 0.01)
{
T = dim(data)[1]
nvar = dim(data)[2]
ydata = data[(lags+1):T,]
xdata = matrix(1, nrow = T - lags, ncol = 1)
for(i in 1:lags){ xdata = cbind(xdata, data[((1+lags)-i):(T-i), ])}
if(!is.null(exogen)) xdata = cbind(xdata, tail(exogen, T - lags))
datamlts = mlts(as.matrix(xdata), as.matrix(ydata), gamma = alpha, ns = ns, nc = nc, delta = delta)
U = datamlts$sigmaR
logdetSigma = log(det(U))
cdelta = (1-delta)/pchisq(qchisq(1 - delta, nvar),nvar+2)
m = sum(datamlts$nooutlier)
phi = nvar*(lags*nvar+1)
ll = -(T-lags)/2 *(log(2*pi)*nvar+log(det(U)))- (m-nvar)*nvar/(2*cdelta)
res = datamlts
res$AIC =-2/(T-lags)*(ll-phi)
res$HQ =-2/(T-lags)*(ll-log(log(T-lags))*phi)
res$SC =-2/(T-lags)*(ll-log(T-lags)*phi/2)
return( res )
}
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.