[R-SIG-Finance] rcoplua.gauss sim problem

[Joe W. Byers]

I am trying to replicate the copulae simulation from the SAS docs using 
gamlss, fCopulae, and copula

My problem is when I try and run rcopula.gauss I get an error message, 
sigma needs to be a correlation matrix in the test for sum(diag(sigma))=d.

No matter what I give rcopula.gauss, correlation matrix or cov matrix I 
get this error.  Even when I pass d=sum(diag(cop.fit$cor)) I get the 
error.  If I comment out the test for a corr matrix the simulation seems 
to work.

Any help is appreciated.

Thank you
Joe


require(VGAM)
require(QRMlib)
require(gamlss)
require(fCopulae)
require(copula)
require(date)

df<-7.5
sig1<-.5
var2<-2.5
seed<-1234
set.seed(seed);
#base on sas code from web site
#generate t distributed data
tdata<-data.frame(date=seq(as.numeric(as.date('06/01/2001')),as.numeric(as.date('11/01/2002')),1))
tdata$tdata<-.05*tdata$date+rTF(length(tdata$date),mu=5000,sigma=sig1,nu=df)

#estimate model data~a * date + constant
     fam<-TF
 
y1<-gamlss(tdata~date+1,data=tdata,family=fam,na.rm=T,nu.start=16.58,sigma.formula=~1
       ,nu.formula=~1)
     y1<-gamlss(tdata~date+1,data=tdata,family=TF(mu.link = 
identity,sigma.link = log, nu.link = log))
     plot(y1,na.rm=T)


  #generate normal garch distributed data
  set.seed(12345)
 
ndata<-data.frame(date=seq(as.numeric(as.date('06/01/2001')),as.numeric(as.date('11/01/2002')),1))
   ndata$rnd<-rnorm(length(ndata$date))
   ndata$v<-NaN
   ndata$e<-NaN
   ndata$v[1]<-.0001+.2*var2^2 +.75* var2
   ndata$e[1]<-sqrt(ndata$v[1])*ndata$rnd[1]
   for (i in 2:length(ndata$date) ) {
     ndata$v[i]<- 0.0001 + 0.2 * ndata$e[(i - 1)]^2 + 0.75 * ndata$v[(i 
- 1)]
     ndata$e[i]<-sqrt(ndata$v[i])+ndata$rnd[i]
   }
   # add constant 25
   ndata$ndata<-25+ndata$e
   # estimate garch 1,1 model of normally distributed data
   yn<-garchFit(ndata~garch(1,1),data=ndata)

  #y = rcauchy(519, loc=-4, scale=1)
  # generate cauchy data
  set.seed(seed);
 
cdata<-data.frame(date=seq(as.numeric(as.date('06/01/2001')),as.numeric(as.date('11/01/2002')),1))
  #cdata$cdata<-unlist(-4+tan(runif(length(cdata$date))-.5)*pi)
  cdata$cdata<-rcauchy(length(cdata$date),loc=-4,scale=1)
  #estimate cauchy model
  yc<-vglm(cdata ~ 1 , data=cdata,cauchy1(lloc="identity"), trace=TRUE, 
crit="c")



  #create transformed uniform residuals for all estimated series' residuals
 
resids<-data.frame(date=tdata$date,tresids=y1$residuals,nresids=yn at residuals,cresids=yc at residuals)
 
resids$tresids<-pnorm(pTF(resids$tresids/exp(y1$sigma.coefficient),y1$nu.coefficients))
  resids$cresids<-pnorm(pcauchy(resids$cresids,yc at coefficients[1]))
 
resids$nresids<-pnorm(resids$nresids,mean(resids$nresids),sd=sd(resids$nresids))
 
yall<-merge(tdata,merge(cdata,ndata,by.x="date",by.y="date",all=T),by.x="date",by.y="date",all=T)
 
rnds<-data.frame(rmvnorm(2000,sigma=cov(resids[-1]),mean=mean(resids[-1])))
  names(rnds)<-c('T',"Normal",'Cauchy')
#fit a normal copula to the uniform residuals
  cop.fit<-fit.norm(resids[,2:4])
#simulate the normal copula
# cop.sim1<-rmnorm(2000, Sigma=cop.fit$Sigma, mu=cop.fit$mu)


#************Here is my problem
  cop.sim<-rcopula.gauss(2000, Sigma=cop.fit$Sigma)

#*******




  #,d=sum(diag(cop.fit$cor)))#, mu=cop.fit$mu)
#transform the simulated correlated normal rv' to uncorrelated 
individual distributions
  sims<-data.frame(cop.sim)
  names(sims)<- c('T',"Normal",'Cauchy')
 
sims$T<-qTF(pnorm(sims$T),mu=y1$mu.coefficients[1],sigma=exp(y1$sigma.coefficients), 
nu=y1$nu.coefficients)
  sims$Normal<-qnorm(sims$Normal,mean(resids$nresids),sd=sd(resids$nresids))
  sims$Cauchy<-qcauchy(pnorm(sims$Cauchy),location=yc at coefficients[1])
 
sims$T<-seq(as.numeric(as.date('11/01/2002'))+1,as.numeric(as.date('11/01/2002'))+2000,1) 
*
    y1$mu.coefficients[2] + sims$T


# tt <- seq(0,10, len=21)                     cauchy1(lloc="loge")
#ncp <- seq(0,6, len=31)
#ptn <- outer(tt,ncp, function(t,d) pt(t, df = 3, ncp=d))
#image(tt,ncp,ptn, zlim=c(0,1),main=t.tit <- "Non-central t - 
Probabilities")
#persp(tt,ncp,ptn, zlim=0:1, r=2, phi=20, theta=200, main=t.tit,
#      xlab = "t", ylab = "non-centrality parameter", zlab = "Pr(T <= t)")


##rcopula.t<-function (n, df, Sigma = equicorr(d, rho), d = 2, rho = 0.7)
##{
##    d <- dim(Sigma)[1]
##    if (sum(diag(Sigma)) != d)
##        stop("Sigma should be correlation matrix")
##    tmp <- rmt(n, df, Sigma)
##    matrix(pt(tmp, df), ncol = d)
##}


#************ rcopula.gauss that works
rcopula.gauss<-function (n, Sigma = equicorr(d, rho), d = 2, rho = 0.7)
{
     d <- dim(Sigma)[1]
#    if (sum(diag(Sigma)) != d)
#        stop("Sigma should be correlation matrix")
     mnorm <- rmnorm(n, Sigma = Sigma)
     matrix(pnorm(mnorm), ncol = d)
}