[R-SIG-Finance] rugarch + VineCopula for value at risk
alexios ghalanos
alexios at 4dscape.com
Fri Aug 22 04:39:42 CEST 2014
Ole,
See a 'corrected' extract of the relevant code section below. You did
not define 'returns_crisis' but I assume this is the same as 'returns'.
You also did not tell us the required packages (please do so in future).
There is absolutely nothing wrong (at least that I can see with the
approach and results).
1. You are filtering the returns by an AR-GARCH model.
2. You estimate the vine copula on the uniform margins based on the
standardized residuals transformed by the conditional distribution from (1).
3. You simulate from this unconditional distribution 100x20 uniform
variates which you transform back to standardized values using the shape
parameter from (1).
3. You then simulate the T+1 conditional distribution (with 100 samples
for each asset) using the GARCH model from (1).
4. You obtain the T+1 distribution of the quantiles for each asset.
Are you very surprised that a particular point in time distribution of
the 5% quantile of certain assets is not entirely in the negative
domain? Think about what it is you are doing, the assumptions, high
degree of uncertainty (you are using 250 points to estimate a GARCH
model) etc The surprise would be if the MEAN of those value was positive.
Am not too familiar with Vine Copulas and the particular package (or if
there is something particular about the way the simulation is done) so I
may have missed something...
-Alexios
Note also the following:
a. For T+1 there is not uncertainty about sigma (see
http://unstarched.net/r-examples/rugarch/a-note-on-simulation-in-the-rugarch-package/
if you are confused about this point).
b. All distributions in rugarch (and their functions qdist, rdist,
pdist, ddist) are based on the standardized (0,1) representation (all
have the location-scale invariance property).
################################################################################
for(i in 1:10)
{
print(i)
windows <- window(returns_crisis[,1:20], start=times[376-250-24+i,1],
end=times[376-25+i,1])
fit <- lapply(windows, ugarchfit, spec=model, solver="hybrid")
print("rugarch fitting done")
residuals <- sapply(fit, residuals, standardize=TRUE)
shape <- sapply(fit, coef)
shape <- shape[5,]
# # I don't know where you found pged, but the correct version is
# pdist("ged",....)
UniformResiduals <- sapply(1:ncol(residuals), function(i){pdist("ged",
residuals[,i], mu=0, sigma=1, shape = shape[i])})
rvine <- RVineStructureSelect(UniformResiduals, indeptest=TRUE,
familyset=familyset)
print(paste(i,"RVine fitting done"))
set.seed(100)
sim <- RVineSim(100, rvine)
print(paste(i,"RVine simulation done"))
for(k in 1:20)
{
# I don't know where you found qged, but the correct version is
# qdist("ged",....)
residuals2 <- qdist("ged", mu=0, sigma=1, sim[,k], shape = shape[k])
# Comment: custom.dist$name must be a character (anything but not NA or
# NULL).
# Also use startMethod="sample" to use the last values of the estimated
# object for the simulation rather than the unconditional values.
rvine_sim <- ugarchsim(fit[[k]], n.sim=1, m.sim=100, custom.dist =
list(name="sample", distfit=matrix(residuals2, ncol=100)),
startMethod="sample")
# cleaner to calculate at each step
rvine_fitted[i,,k] <- as.numeric(fitted(rvine_sim))
rvine_sigma[i,,k] <- as.numeric(sigma(rvine_sim))
VaR01[i,,k] <- rvine_fitted[i,,k] + rvine_sigma[i,,k] * qdist('ged',
0.01, mu=0, sigma=1, shape = shape[k])
VaR05[i,,k] <- rvine_fitted[i,,k] + rvine_sigma[i,,k] * qdist('ged',
0.05, mu=0, sigma=1, shape = shape[k])
}
print(paste(i, "rugarch forecast done"))
}
################################################################################
On 22/08/2014 00:44, Ole Bueker wrote:
> Sorry for sending another message, but it seems like there is still an
> error in the formatting of the code..
> I'll share the R code using dropbox, in case the code is still not
> readable in this
> email: https://www.dropbox.com/s/rc9kpanule52zu8/minimum_working.R
>
> #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_crisis))
> windows <- matrix(0, 20, 250)
> familyset <- c(1:5, 7, 10, 13, 14, 17, 20)
> sim <- array(0, dim = c(100, 20))
> residuals2 <- array(0, dim = c(100, 20))
> rvine_fitted <- array(0, dim = c(10,100,20))
> rvine_sigma <- array(0, dim = c(10,100,20))
> VaR01 = VaR05 = array(0, dim = c(10,100,20))
>
> #Main calculation
> for(i in 1:10)
> {
> print(i)
> windows <- window(returns_crisis[,1:20], start=times[376-250-24+i,1],
> end=times[376-25+i,1]) # Step 1
> fit <- lapply(windows, ugarchfit, spec=model, solver="hybrid")
> print("rugarch fitting done")
>
> residuals <- sapply(fit, residuals, standardize=TRUE) # Step 2
> shape <- sapply(fit, coef)
> shape <- shape[5,]
>
> UniformResiduals <- pged(residuals, nu = shape) # Step 3
> 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) # Step 4
> print(paste(i,"RVine fitting done"))
> for(j in 1:100)
> {
> sim[j,] <- RVineSim(1, rvine) # Step 5
> }
> print(paste(i,"RVine simulation done"))
> for(k in 1:20)
> {
> residuals2[,] <- qged(sim[,], nu = shape[k]) # Step 6
>
> residuals_temp <- residuals2[,k]
> rvine_sim <- ugarchsim(fit[[k]], n.sim=1, m.sim=100, custom.dist =
> list(name=NA, distfit=residuals_temp)) # Step 7
>
> rvine_fitted[i,,k] <- fitted(rvine_sim) # Step 8
> rvine_sigma[i,,k] <- sigma(rvine_sim)
> if (i==10)
> {
> for(j in 1:100)
> {
> VaR01[,j,k] <- rvine_fitted[,j,k] + rvine_sigma[,j,k] * qdist('ged',
> 0.01, mu=0, sigma=1, shape = shape[k]) # Step 9
> VaR05[,j,k] <- rvine_fitted[,j,k] + rvine_sigma[,j,k] * qdist('ged',
> 0.05, mu=0, sigma=1, shape = shape[k])
> }
> }else {}
> }
> print(paste(i, "rugarch forecast done"))
> }
>
>
> #Cleanup
> remove(i, j, k, ix, familyset, model, sim, rvine_sim, residuals,
> residuals_temp)
>
>
> ------------------------------------------------------------------------
>
>> From: ole.bueker at outlook.com
>> To: alexios at 4dscape.com; r-sig-finance at r-project.org
>> Date: Fri, 22 Aug 2014 01:10:55 +0200
>> Subject: Re: [R-SIG-Finance] rugarch + VineCopula for value at risk
>>
>>
>>
>>
>> I have no idea what happened to the formatting, looked fine to me when
> I typed it in Outlook.
>> Here's the revised version, with amended code so that it is minimally
> reproducible (for this purpose, I will only use the first 25 companies,
> 100 simulations for 10 trading days).
>> The revised code should take around 5 minutes,
>>
>>
>> Again, heres the code summary:
>> 1. Fit GARCH models to each series.
>> 2. Extract standardized returns.
>> 3. Transform standardized returns to uniform marginals using the
> parametric IFM method by Joe.
>> 4. Fit vine copulas.
>> 5. Generate 100 1-day ahead forecasts from the vine copulas.
>> 6. Reverse transform the simulated values.
>> 7. Use these transformed forecasts in ugarchsim (using custom.dist)
>> 8. Extract forecasted mu and sigma.
>> 9. Calculate 95% and 99% VaR.
>>
>> The full code:
>> #Load data and define variablesreturns <- 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_crisis))windows <- matrix(0, 20,
> 250)familyset <- c(1:5, 7, 10, 13, 14, 17, 20)sim <- array(0, dim =
> c(100, 20))residuals2 <- array(0, dim = c(100, 20))rvine_fitted <-
> array(0, dim = c(10,100,20))rvine_sigma <- array(0, dim =
> c(10,100,20))VaR01 = VaR05 = array(0, dim = c(10,100,20))
>> #Main calculation
>> for(i in 1:10){ print(i) windows <- window(returns_crisis[,1:20],
> start=times[376-250-24+i,1], end=times[376-25+i,1]) # Step 1 fit <-
> lapply(windows, ugarchfit, spec=model, solver="hybrid") print("rugarch
> fitting done")
>> residuals <- sapply(fit, residuals, standardize=TRUE) # Step 2 shape
> <- sapply(fit, coef) shape <- shape[5,]
>> UniformResiduals <- pged(residuals, nu = shape) # Step 3
> 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) # Step 4 print(paste(i,"RVine fitting done")) for(j
> in 1:100) { sim[j,] <- RVineSim(1, rvine) # Step 5 }
> print(paste(i,"RVine simulation done")) for(k in 1:20) { residuals2[,]
> <- qged(sim[,], nu = shape[k]) # Step 6
>> residuals_temp <- residuals2[,k] rvine_sim <- ugarchsim(fit[[k]],
> n.sim=1, m.sim=100, custom.dist = list(name=NA, distfit=residuals_temp))
> # Step 7
>> rvine_fitted[i,,k] <- fitted(rvine_sim) # Step 8 rvine_sigma[i,,k] <-
> sigma(rvine_sim) if (i==10) { for(j in 1:100) { VaR01[,j,k] <-
> rvine_fitted[,j,k] + rvine_sigma[,j,k] * qdist('ged', 0.01, mu=0,
> sigma=1, shape = shape[k]) # Step 9 VaR05[,j,k] <- rvine_fitted[,j,k] +
> rvine_sigma[,j,k] * qdist('ged', 0.05, mu=0, sigma=1, shape = shape[k])
> } }else {} } print(paste(i, "rugarch forecast done"))}
>>
>> #Cleanupremove(i, j, k, ix, familyset, model, sim, rvine_sim,
> residuals, residuals_temp)
>>
>> Hope this time everything is formatted correctly!
>> After running the code, I export VaR01 and VaR05 to Excel, and notice
> that there's usually quite a few positive values.
>> This seems to happen because sometimes the sigma is too low compare to
> the (positive) forecast, and therefore the VaR stays in the positive.
>> My guess is that either there's a mistake in my code at step 3 (maybe
> the wrong shape is used?).Or more likely there's a mistake at step 6/7
> (wrong transformation or transformed values are not inserted correctly
> into ugarchsim).
>> > Date: Thu, 21 Aug 2014 18:21:49 +0100
>> > From: alexios at 4dscape.com
>> > To: ole.bueker at outlook.com; r-sig-finance at r-project.org
>> > Subject: Re: [R-SIG-Finance] rugarch + VineCopula for value at risk
>> >
>> > Ole,
>> >
>> > The email you've sent is really badly formatted and all over the place.
>> >
>> > 1. Try to send a well structured text-only email (with no special
>> > characters).
>> > 2. Make an effort to send a minimally reproducible example.
>> > Telling us that your problem took 8 hours to run is not the way to go.
>> > Take a small subsample of your data which generates the behavior you
>> > want to highlight (try to find one that does), and use that so that we
>> > can investigate with minimal fuss.
>> >
>> > -Alexios
>> >
>> > On 21/08/2014 17:50, Ole Bueker wrote:
>> > > 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 si
>> > mulations.
>> > > 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]]
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-SIG-Finance at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R
> questions should go.
More information about the R-SIG-Finance
mailing list