[R-SIG-Finance] rcopula.gauss problem with sum(diag(correlationmatrix))

[Joe W. Byers]

>
>

All,

I keep getting rcopula.gauss to error out due to machine/R significant digits. 
The test in the methods sum(diag(sigma))==d fails because the one of my diag
elements is not quite 1.  I had to execute a ceiling function on the diag
elements of the correlation matrix to get the method not to fail.  

Is this a problem for anyone else?

My example is
              
              
              require(VGAM)
              require(QRMlib)
              require(gamlss)
              require(fCopulae)
              require(fGarch)
              require(copula)
              require(date)
              require(rgl)
              
              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+.25*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]
                }
                gs<-garchSim(model = list( mu=25,omega = .0001, alpha = .2, 
                beta = .75), n = length(ndata$date), 
                  n.start = 100, presample = NULL, cond.dist = c("rnorm"), 
                  rseed = 12345)
              
                # add constant 25
                ndata$ndata<-as.numeric(gs); #25+ndata$e
                # estimate garch 1,1 model of normally distributed data
                yn<-garchFit(ndata~garch(1,1),data=ndata,trace=F)
               
               #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))
               cauch<-unlist(-4+tan(runif(length(cdata$date))-.5)*pi)
               set.seed(seed);
               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)
               
 ## to get rcopula.gauss to not fail
# diag(cop.fit$cor)=ceil(diag(cop.fit$cor))
 #uncomment the above line to get rcopula.gauss to work
 cop.sim<-rcopula.gauss(2000, Sigma=cop.fit$cor)
 
           ####################
           sims<-data.frame(cop.sim)
           names(sims)<- c('T',"Normal",'Cauchy')
           # for rmnorm
           #sims$T<-qTF(pnorm(sims$T),mu=y1$mu.coefficients[1],
           # sigma=exp(y1$sigma.coefficients), nu=y1$nu.coefficients)
           #sims$Normal<-qnorm(pnorm(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
           # for rcopula.gauss
           sims$T<-qTF(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(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
          
          f2<-kde2d(sims$T,sims$Cauchy,lims=c(range(sims$T),c(-150,150)))
          open3d()
          bg3d("white")
          material3d(col="black")
           persp3d(x=f2$x,y=f2$y,z=f2$z, aspect=c(1, 1, 0.5), col = "lightblue",
            xlab='T',ylab='Cauchy')