Problems in computing VaR with GARCH-GPD-copula approach

Question

I use a time-varying Gaussian copula (with GARCH-filtered standardized residuals modeled semiparametrically with Gaussian kernel interior and GPD tails, i.e. generalized pareto distributed) to simulate N copula realizations and convert them back into portfolio returns (equally-weighted). Using a testing windows of 2000 daily returns, I get a number of VaR violations which exceeds the expected number of violations by far (e.g. 450 instead of 200 for 90%-VaR, 280 instead of 50 for 97.5-VaR, ...). However, if I analyze the simulated returns of each component individually, each stock has a number of violations which is approximately in line with those expected. Since the GARCH-GPD-copula approach is quite established in the literature for risk management purposes, my results are most probably due to some mistake I've made, I exclude (for the moment, at least) that I have opposing findings with respect to those usually published.

The only two options that come to my mind are: i) effect of N; ii) code doesn't work properly (due to my inexperience with R).

i) I've simulated 5000 correlated copula realizations. Is 5000 maybe not enough, would accuracy maybe increase if N=10'000? Is there any rule about how to select number of simulations?

ii) This is the essential part of the code (T equals 2000)

tmp=cgarchfilter(specx, returns[1:2000,], filter.control=list(n.old=T))
presigma=matrix(tail(sigma(tmp),1),ncol=9)
prereturns=matrix(unlist(returns[T+i-1,]),ncol=9,nrow=1)
preresiduals=matrix(tail(residuals(tmp),1),ncol=9, nrow=1)
preR=last(rcor(tmp))[,,1]
diag(preR)=1
preQ=tmp@mfilter$Qt[[length(tmp@mfilter$Qt)]]
preZ=tail(tmp@mfilter$Z,1)

for (i in 1:2000){
tmp=cgarchfilter(specx, returns[i:(T+i-1),], filter.control=list(n.old=T))
sim <- cgarchsim(fit, n.sim=1, m.sim=5000, startMethod="sample", preR=preR, preQ=preQ,
               preZ=preZ, prereturns=prereturns, preresiduals=preresiduals, presigma=presigma,
               cluster=NULL)
simx=t(sapply(sim@msim$simX, FUN=function(x) x[1,]))
a[i,,] <- simx
}

specx fixates the GARCH-GPD parameters of the fitting, all pre_ objects are needed as inputs in the cgarchsim function. All inputs outside the loop compute the inputs for i=1; these are then reevaluated in the loop from i=1 to i=2000 (in-sample and testing periods are both over 2000 data, 4000 in total). simx then contains all 5000 1-step ahead simulated returns. So a is a 2000x9x5000 array (9 is the number of series).

Do you know what could have gone wrong? I'm lost.

EDIT 1: It has to be the code. I've tried with N=10'000 and the number of violations doesn't change a bit. I've also tried considering smaller return intervals (e.g. last 100 instead of last 2000) to compute all pre_ objects and account more heavily for recent news, and still the number of violations is the same. Am I wrong (which I assume), or am I swimming in an undiscovered (and flawed) area of the GARCH-copula sea?