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@residuals,cresids=yc@residuals)
resids$tresids<-pnorm(pTF(resids$tresids/
exp(y1$sigma.coefficient),y1$nu.coefficients))
resids$cresids<-pnorm(pcauchy(resids$cresids,yc@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@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@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')
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@residuals,cresids=yc@residuals)
resids$tresids<-pnorm(pTF(resids$tresids/
exp(y1$sigma.coefficient),y1$nu.coefficients))
resids$cresids<-pnorm(pcauchy(resids$cresids,yc@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@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@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')
No comments:
Post a Comment
Thank you