[R-SIG-Finance] rugarch + VineCopula for value at risk

[Ole Bueker]

Hello,
I am trying to use the R packages rugarch and VineCopula for simulating returns of 112 companies for a time period of 25 days with daily re-estimations. After the simulation, I wish to calculate the 99% and 95% value-at-risk and compare them to the actual returns.I use a moving window of 250 days and 1000 simulations per iteration.
(This estimation is quite time-intensive as 112 companies might be too many for VineCopula..)
The overall loop takes around 8 hours on my home computer, so I wouldn’t recommend to actually run the code. My problem is that some of my calculated value at risk forecasts seem to be positive – this is not a “just invert the VaR” kind of problem (at least I don’t think so). 
To summarize the code:
1.  Fit GARCH models to each series.2.  Extract standardized returns (and shape parameters)3. Transform standardized returns to uniform marginals using parametric method (IFM by Joe, 1987).4. Fit vine copulas5. Generate a 1000 x 112  matrix (1000 1-day ahead forecasts for all 112 companies)6. Reverse transform the simulated values.7. Use these transformed forecasts in ugarchsim
8. Extract forecasted values & sigmas.9. Calculate Value-at-Risk. 
Anyway, here’s my code so far:# Load Data and define variables
returns <-  read.zoo("E:/Dropbox/my own/Programming/R/returns.csv", header=TRUE, sep=",", format="%d-%m-%y")model<-ugarchspec(variance.model=list(model="sGARCH",garchOrder=c(1,1)),mean.model=list(armaOrder=c(1,0),include.mean=FALSE),distribution.model="ged")times <- as.data.frame(time(returns))windows <- matrix(0, 112, 250)familyset <- c(1:5, 7, 10, 13, 14, 17, 20)                    # The vine copulas to be testedsim <- array(0, dim = c(1000, 112))residuals2 <- array(0, dim = c(1000, 112))rvine_fitted <- array(0, dim = c(25,1000,112))rvine_sigma <- array(0, dim = c(25,1000,112))VaR01 = VaR05 = array(0, dim = c(25,1000,112))

#Main calculation
for(i in 1:25){  print(i)  windows <- window(returns_crisis, start=times[376-250-24+i,1], end=times[376-25+i,1])            #Define the moving window  fit <- lapply(windows, ugarchfit, spec=model, solver="hybrid")                                                                #Fit the garch models  print("rugarch fitting done")  residuals <- sapply(fit, residuals, standardize=TRUE)                                                                                     #Extract residuals & shape parameters  shape <- sapply(fit, coef)  shape <- shape[5,]  UniformResiduals <- pged(residuals, nu = shape)                                                                                           #Transform residuals into uniform marginals  if(any(UniformResiduals > 0.99999))  {    ix = which(UniformResiduals > 0.99999)   UniformResiduals [ix] = 0.99999  }  if(any(UniformResiduals < .Machine$double.eps))  {    ix = which(UniformResiduals < (1.5*.Machine$double.eps))    UniformResiduals [ix] = .Machine$double.eps  }  rvine <- RVineStructureSelect(UniformResiduals, indeptest=TRUE, familyset=familyset)             #Fit the Vine copulas  print(paste(i,"RVine fitting done"))  for(j in 1:1000)                                                                                                                                                                 #Simulate 1000 1-day ahead using VineCopula  {  sim[j,] <- RVineSim(1, rvine)                                                                                                                                     # 1000 x 112 matrix of forecasts  }print(paste(i,"RVine simulation done"))                 for(k in 1:112)                                                                                                                                                    #Next: ugarchsimulation for all 112 companies                {                residuals2[,] <- qged(sim[,], nu = shape[k])                                                                                         # 1000 x 112 matrix of standardized residuals                residuals_temp <- residuals2[,k]                                                                                                               # 1000 x 1 vector of standardized residuals for individual company                rvine_sim <- ugarchsim(fit[[k]], n.sim=1, m.sim=1000, custom.dist = list(name=NA, distfit=residuals_temp))  #1000 simulations using the standardized residuals from Vine copula models for ugarchfit                rvine_fitted[i,,k] <- fitted(rvine_sim)                                                                                                     #Extract forecasted values - 25 x 1000 x 112               rvine_sigma[i,,k] <- sigma(rvine_sim)                                                                                                     #Extract forecasted sigmas - 25 x 1000 x 112               for(j in 1:1000)                                                                                                                                                                   #Next: Value at risk              {                VaR01[,j,k] <- rvine_fitted[,j,k] + rvine_sigma[,j,k] * qdist('ged', 0.01, mu=0, sigma=1, shape = shape[k]) #Value at risk for 99% quantile                VaR05[,j,k] <- rvine_fitted[,j,k] + rvine_sigma[,j,k] * qdist('ged', 0.05, mu=0, sigma=1, shape = shape[k]) #Value at risk for 95% quantile               }                 }}remove(i, j, k)                                                                                                                                                                   #Cleanupremove(windows, fit, residuals, shape, residuals2, residuals_temp, rvine, sim, rvine_sim)            #Cleanup   Hope I didn’t make any mistakes in my approach, but it seems like this is the “standard” copula + rugarch approach – if anyone is familiar with this, I am open to suggestions on how to speed up the simulations.
So far, so good – the problem I am facing now:
Some (only a few) of my value at risk values are positive..I have manually checked and it seems like the fitted value is much larger than the sigma, so Value at Risk is positive – which doesn’t really make any economic sense to me.
Here’s a dropbox link to the returns.csv, in case anyone is interested in running my code: https://www.dropbox.com/s/69i5959f3h4kweb/returns.csv

Best Regards,
Ole 		 	   		  
	[[alternative HTML version deleted]]