R/JGQD.density.R
In eta21/DiffusionRjgqd: Inference and Analysis for Jump Generalized Quadratic Diffusions
globalVariables(c('mua','alph'))
JGQD.density=function(Xs=4,Xt=seq(5,8,1/10),s=0,t=5,delt=1/100,Jdist='Normal',Jtype='Add',Dtype='Saddlepoint',Trunc=c(8,4),factorize=FALSE,factor.type='Diffusion',beta,print.output=TRUE,eval.density=TRUE)
{
Integ = TRUE
P=100;alpha=0.1;lower=0;upper=50
TR.order=Trunc[1]
DTR.order=Trunc[2]
Dtypes =c('Saddlepoint','Edgeworth','Normal.A','Gamma.A','Beta','Saddle.Gamma','Gamma.Saddle')
check_for_model=function()
{
txt=''
namess=c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta','Ja','Jb')
func.list=rep(0,length(namess))
obs=objects(pos=1)
for(i in 1:length(namess))
{
if(sum(obs==namess[i])){func.list[i]=1}
}
check=F
if(sum(func.list)==0)
{
txt='
--------------------------------------------------------------------------------
No model has been defined yet! Try for example:
--------------------------------------------------------------------------------
GQD.remove()
G0=function(t){5+2*sin(4*pi*t)}
G1=function(t){-1}
Q1=function(t){1}
Lam0 = function(t){0.1}
Jmu = function(t){0}
Jsig = function(t){1}
model=JGQD.density(5,seq(0,10,1/10),0,5,1/100)
JGQD.plot(model)
--------------------------------------------------------------------------------
'
check=T
}
if((sum(func.list)>0)&&(sum(func.list[c(1:3)+3])==0))
{
txt2='
--------------------------------------------------------------------------------
No diffusion coefficients specified. The resulting process will be that of a
Hawkes process.
{factor.type} changed to "Hawkes"
--------------------------------------------------------------------------------
'
warning(txt2)
factor.type = 'Hawkes'
}
return(list(check=check,txt=txt))
}
check_for=check_for_model()
if(check_for[[1]]){stop(check_for[[2]])}
TR.order=Trunc[1]
DTR.order=Trunc[2]
#Dtypes =c('Saddlepoint','Normal','Beta')
Dindex = which(Dtypes==Dtype)
JDtypes =c('Normal','Exponential','Gamma','Laplace','Uniform')
JDindex = which(JDtypes==Jdist)
b1 = '\n==============================================================================\n'
b2 = '==============================================================================\n'
warn=c(
'Input 1: Argument {lower} must be of length 1.\n'
,'Input 2: Argument {upper} must be of length 1.\n'
,'Input 3: {upper} must be > {lower}.\n'
,'Input 4: Range [lower,upper] must be strictly positive for Dtype Gamma or InvGamma.\n'
,'Input 5: Dtype cannot be Beta for observations not in (0,1).\n'
,'Input 6: P must be >= 10.\n'
,'Input 7: Dtype has to be one of Saddlepoint or Edgeworth.\n'
,'Input 8: Trunc[2] must be <= Trunc[1].\n'
,'Input 9: Large {delt} may result in poor approximations.\n'
,'Input 10: Argument {delt} must be < 1.\n'
,'Input 11: length(P)!=1.\n'
,'Input 12: length(delt)!=1.\n'
,'Input 13: length(alpha)!=1.\n'
,'Input 14: length(Xt) must be > 1.\n'
,'Input 15: Argument {Xs} must be of length 1.\n'
,'Input 16: Argument {Xt} must be of type vector!.\n'
,'Input 17: Argument {time} must be of type vector!.\n'
,'Input 18: Argument {alpha} must be of length 1.\n'
,'Input 19: Argument {Trunc} must be of length 2.\n'
,'Input 20: Cumulant truncation (Trunc[1]) must be either 4 or 8.\n'
,'Input 21: Density truncation (Trunc[2]) must be either 4 or 8.\n'
,'Input 22: Starting time {s} cannot be greater than {t}.\n'
,'Input 23: {Jdist} must be one of Normal, Exponential, Gamma or Laplace.\n'
,'Input 24: {Jtype} has to be of type Add or Mult.\n'
,'Input 25: {factorize} has to be TRUE or FALSE.\n'
,'Input 26: {factor.type} has to be Diffusion or Hawkes.\n'
,'Input 27: Diffusion terms are 0. Use {factorize=TRUE} and {factor.type="Hawkes"}.\n'
,'Input 28: Diffusion terms are 0. Density approximant cab perform poorly without \nfactorization. Try {factorize=TRUE}.\n'
)
warntrue=rep(FALSE,40)
if(length(Xt)<2){warntrue[14] =TRUE}#{stop(paste0(b1,warn[14],b2))}
if(length(Xs)!=1){warntrue[15] =TRUE}#{stop(paste0(b1,warn[15],b2))}
if(!is.vector(Xt)){warntrue[16] =TRUE}#{stop(paste0(b1,warn[16],b2))}
if(sum(Dindex)==0){warntrue[7] =TRUE}#{stop(paste0(b1,warn[7],b2))}
if(length(lower)>1){warntrue[1] =TRUE}#{stop(paste0(b1,warn[1],b2))}
if(length(upper)>1){warntrue[2] =TRUE}#{stop(paste0(b1,warn[2],b2))}
if(upper<=lower) {warntrue[3] =TRUE}# {stop(paste0(b1,warn[3],b2))}
#if((Dindex==3)|(Dindex==4)){if(lower[1]<=0){warntrue[4] =TRUE}}#{stop(paste0(b1,warn[4],b2))}}
if(Dindex==5){if(any(Xt<=0)|any(Xt>=1)){warntrue[5] =TRUE}}#{stop(paste0(b1,warn[5],b2))}}
if(P<10){warntrue[6] =TRUE}#{stop(paste0(b1,warn[6],b2))}
if(length(P)!=1){warntrue[11] =TRUE}#{stop(paste0(b1,warn[11],b2))}
if(length(alpha)!=1){warntrue[18] =TRUE}#{stop(paste0(b1,warn[18],b2))}
if(Trunc[2]>Trunc[1]){warntrue[8] =TRUE}#{stop(paste0(b1,warn[8],b2))}
if(length(Trunc)!=2){warntrue[19] =TRUE}#{stop(paste0(b1,warn[19],b2))}
if(sum(TR.order==c(4,8))!=1){warntrue[20] =TRUE}#{stop(paste0(b1,warn[20],b2))}
if(sum(DTR.order==c(4,8))!=1){warntrue[21] =TRUE}#{stop(paste0(b1,warn[21],b2))}
#if(delt>=1/5){stop(paste0(b1,warn[9],b2))}
if(delt>=1){warntrue[10] =TRUE}#{stop(paste0(b1,warn[10],b2))}
if(t<s){warntrue[22] =TRUE}#{stop(paste0(b1,warn[22],b2))}
if((Jtype!='Add')&&(Jtype!='Mult')) {warntrue[24] =TRUE}
if((factorize!=TRUE)&&(factorize!=FALSE)) {warntrue[25] =TRUE}
if((factor.type!='Diffusion')&&(factor.type!='Hawkes')) {warntrue[26] =TRUE}
if(sum(JDindex)==0){stop(paste0(b1,warn[23],b2))}
# Print output:
if(any(warntrue))
{
prnt = b1
for(i in which(warntrue))
{
prnt = paste0(prnt,warn[i])
}
prnt = paste0(prnt,b2)
stop(prnt)
}
diffusion.function=function()
{
time =seq(s,t,delt)[-1]
bod='0*x^0'
coefs=c('Q0','Q1','Q2')
obs=objects(pos=1)
for(i in 1:length(coefs))
{
if(sum(obs==coefs[i])){bod=paste0(bod,'+(',body(coefs[i])[2],')*x^',i-1)}
}
difff=function(x,t){}
body(difff) = parse(text=bod)
#print(body(difff))
M = matrix(0,length(Xt),length(time))
for(i in 1:length(time))
{
M[,i]=round(difff(Xt,time[i]),6)>0
}
if(sum(M)==length(time)*length(Xt)){return(1)}
if(sum(M)!=0) {return(2)}
if(sum(M)==0) {return(3)}
}
switch(diffusion.function(),
{},
{},
{
if((factorize)&&(factor.type!='Hawkes')){stop(paste0(b1,warn[27],b2))}
if(!factorize){warning(paste0(b1,warn[28],b2))}
})
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
function.list=objects(pos=1)
namess=c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta','Ja','Jb')
coef.terms2=paste0(namess,'(t)')
checklist=rep(0,length(namess))
for(i in 1:length(namess))
{
checklist[i]= sum(function.list== namess[i])
}
if(factorize)
{
if(sum(checklist[c(3,5,6)])==0){Dtype='Normal.A'}
}
if(Jdist=='Normal') {checklist[12:16]=0}
if(Jdist=='Exponential'){checklist[c(10,11,13,14,15,16)]=0}
if(Jdist=='Gamma') {checklist[c(10:12,15,16)]=0}
if(Jdist=='Laplace') {checklist[c(10:14)]=0}
if(Jdist=='Uniform') {checklist[c(10:14)]=0}
coef.index=which(checklist==1)
#==============================================================================
# Interface Module
#==============================================================================
namess4=namess
trim <- function (x) gsub("([[:space:]])", "", x)
for(i in 1:length(namess4))
{
if(sum(function.list==namess4[i]))
{
namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
}
}
namess4=matrix(namess4,length(namess4),1)
dinfo = c('Density approx. : ',
'Trunc. Order : ',
'Dens. Order : ')
dinfo[1] =paste0(dinfo[1],Dtype)
#if(Dindex!=1)
#{
# dinfo[2] =paste0(dinfo[2],P)
# dinfo[3] =paste0(dinfo[3],alpha)
#}
dinfo[2] =paste0(dinfo[2],Trunc[1])
dinfo[3] =paste0(dinfo[3],Trunc[2])
buffer0=c('================================================================')
buffer1=c('----------------------------------------------------------------')
buffer2=c('................................................................')
buffer3=c('... ... ... ... ... ... ... ... ... ... ... ')
buffer4=c('_____________________ Drift Coefficients _______________________')
buffer5=c('___________________ Diffusion Coefficients _____________________')
buffer6=c('__________________ Distribution Approximant ____________________')
buffer7=c('_______________________ Jump Mechanism _________________________')
buffer8=c('......................... Intensity ............................')
buffer9=c('........................... Jumps ..............................')
type.sol =" Jump Generalized Quadratic Diffusion (JGQD) "
Info=c(buffer0,type.sol,buffer0,buffer4,namess4[1:3],buffer5,namess4[4:6],buffer7,buffer8,namess4[7:9],buffer9,Jdist,namess4[9+which(checklist[-(1:9)]==1)],buffer6,dinfo)
Info=data.frame(matrix(Info,length(Info),1))
colnames(Info)=''
if(print.output)
{
print(Info,row.names = FALSE,right=F)
}
if(Integ)
{
if(TR.order ==8)
{
MAT=rbind(
c('1','1*x[1]','1*x[2]','','',''),
c('2*x[1]','2*x[2]','2*x[3]','1','1*x[1]','1*x[2]'),
c('3*x[2]','3*x[3]','3*x[4]','3*x[1]','3*x[2]','3*x[3]'),
c('4*x[3]','4*x[4]','4*x[5]','6*x[2]','6*x[3]','6*x[4]'),
c('5*x[4]','5*x[5]','5*x[6]','10*x[3]','10*x[4]','10*x[5]'),
c('6*x[5]','6*x[6]','6*x[7]','15*x[4]','15*x[5]','15*x[6]'),
c('7*x[6]','7*x[7]','7*x[8]','21*x[5]','21*x[6]','21*x[7]'),
c('8*x[7]','8*x[8]','8*(-(-9*x[1]*x[8]-36*x[2]*x[7]+72*x[1]^2*x[7]-84*x[3]*x[6]+504*x[1]*x[2]*x[6]-504*x[1]^3*x[6]-126*x[4]*x[5]+1008*x[1]*x[3]*x[5]+756*
x[2]^2*x[5]-4536*x[1]^2*x[2]*x[5]+3024*x[1]^4*x[5]+630*x[1]*x[4]^2+2520*x[2]*x[3]*x[4]-7560*x[1]^2*x[3]*x[4]-11340*x[1]*x[2]^2*x[4]+30240*x[1]^3*x[2]*x[4]
-15120*x[1]^5*x[4]+560*x[3]^3-15120*x[1]*x[2]*x[3]^2+20160*x[1]^3*x[3]^2-7560*x[2]^3*x[3]+90720*x[1]^2*x[2]^2*x[3]-151200*x[1]^4*x[2]*x[3]+60480*
x[1]^6*x[3]+22680*x[1]*x[2]^4-151200*x[1]^3*x[2]^3+272160*x[1]^5*x[2]^2-181440*x[1]^7*x[2]+40320*x[1]^9))','28*x[6]','28*x[7]','28*x[8]'))
MAT2=rbind(
c('1','1*x[8+1]','1*x[8+2]','','',''),
c('2*x[8+1]','2*x[8+2]','2*x[8+3]','1','1*x[8+1]','1*x[8+2]'),
c('3*x[8+2]','3*x[8+3]','3*x[8+4]','3*x[8+1]','3*x[8+2]','3*x[8+3]'),
c('4*x[8+3]','4*x[8+4]','4*x[8+5]','6*x[8+2]','6*x[8+3]','6*x[8+4]'),
c('5*x[8+4]','5*x[8+5]','5*x[8+6]','10*x[8+3]','10*x[8+4]','10*x[8+5]'),
c('6*x[8+5]','6*x[8+6]','6*x[8+7]','15*x[8+4]','15*x[8+5]','15*x[8+6]'),
c('7*x[8+6]','7*x[8+7]','7*x[8+8]','21*x[8+5]','21*x[8+6]','21*x[8+7]'),
c('8*x[8+7]','8*x[8+8]','8*(-(-9*x[8+1]*x[8+8]-36*x[8+2]*x[8+7]+72*x[8+1]^2*x[8+7]-84*x[8+3]*x[8+6]+504*x[8+1]*x[8+2]*x[8+6]-504*x[8+1]^3*x[8+6]-126*x[8+4]*x[8+5]+1008*x[8+1]*x[8+3]*x[8+5]+756*
x[8+2]^2*x[8+5]-4536*x[8+1]^2*x[8+2]*x[8+5]+3024*x[8+1]^4*x[8+5]+630*x[8+1]*x[8+4]^2+2520*x[8+2]*x[8+3]*x[8+4]-7560*x[8+1]^2*x[8+3]*x[8+4]-11340*x[8+1]*x[8+2]^2*x[8+4]+30240*x[8+1]^3*x[8+2]*x[8+4]
-15120*x[8+1]^5*x[8+4]+560*x[8+3]^3-15120*x[8+1]*x[8+2]*x[8+3]^2+20160*x[8+1]^3*x[8+3]^2-7560*x[8+2]^3*x[8+3]+90720*x[8+1]^2*x[8+2]^2*x[8+3]-151200*x[8+1]^4*x[8+2]*x[8+3]+60480*
x[8+1]^6*x[8+3]+22680*x[8+1]*x[8+2]^4-151200*x[8+1]^3*x[8+2]^3+272160*x[8+1]^5*x[8+2]^2-181440*x[8+1]^7*x[8+2]+40320*x[8+1]^9))','28*x[8+6]','28*x[8+7]','28*x[8+8]'))
dims=rep('',8)
dims2=rep('',8)
for(i in which(checklist[1:6]==1))
{
for(j in 1:8)
{
if(MAT[j,i]!='')
{
dims[j] = paste0(dims[j],'+(',body(namess[i])[2],')*',MAT[j,i])
dims2[j] = paste0(dims2[j],'+(',body(namess[i])[2],')*',MAT2[j,i])
}
}
}
if(Jtype=='Add')
{
JA0 =
c("+1*mm1"
,"+2*mm1*x[1]+1*mm2"
,"+3*mm1*x[2]+3*mm2*x[1]+1*mm3"
,"+4*mm1*x[3]+6*mm2*x[2]+4*mm3*x[1]+1*mm4"
,"+5*mm1*x[4]+10*mm2*x[3]+10*mm3*x[2]+5*mm4*x[1]+1*mm5"
,"+6*mm1*x[5]+15*mm2*x[4]+20*mm3*x[3]+15*mm4*x[2]+6*mm5*x[1]+1*mm6"
,"+7*mm1*x[6]+21*mm2*x[5]+35*mm3*x[4]+35*mm4*x[3]+21*mm5*x[2]+7*mm6*x[1]+1*mm7"
,"+8*mm1*x[7]+28*mm2*x[6]+56*mm3*x[5]+70*mm4*x[4]+56*mm5*x[3]+28*mm6*x[2]+8*mm7*x[1]+1*mm8")
JA1 =
c("+1*mm1*x[1]"
,"+2*mm1*x[2]+1*mm2*x[1]"
,"+3*mm1*x[3]+3*mm2*x[2]+1*mm3*x[1]"
,"+4*mm1*x[4]+6*mm2*x[3]+4*mm3*x[2]+1*mm4*x[1]"
,"+5*mm1*x[5]+10*mm2*x[4]+10*mm3*x[3]+5*mm4*x[2]+1*mm5*x[1]"
,"+6*mm1*x[6]+15*mm2*x[5]+20*mm3*x[4]+15*mm4*x[3]+6*mm5*x[2]+1*mm6*x[1]"
,"+7*mm1*x[7]+21*mm2*x[6]+35*mm3*x[5]+35*mm4*x[4]+21*mm5*x[3]+7*mm6*x[2]+1*mm7*x[1]"
,"+8*mm1*x[8]+28*mm2*x[7]+56*mm3*x[6]+70*mm4*x[5]+56*mm5*x[4]+28*mm6*x[3]+8*mm7*x[2]+1*mm8*x[1]")
JA2 =
c("+1*mm1*x[2]"
,"+2*mm1*x[3]+1*mm2*x[2]"
,"+3*mm1*x[4]+3*mm2*x[3]+1*mm3*x[2]"
,"+4*mm1*x[5]+6*mm2*x[4]+4*mm3*x[3]+1*mm4*x[2]"
,"+5*mm1*x[6]+10*mm2*x[5]+10*mm3*x[4]+5*mm4*x[3]+1*mm5*x[2]"
,"+6*mm1*x[7]+15*mm2*x[6]+20*mm3*x[5]+15*mm4*x[4]+6*mm5*x[3]+1*mm6*x[2]"
,"+7*mm1*x[8]+21*mm2*x[7]+35*mm3*x[6]+35*mm4*x[5]+21*mm5*x[4]+7*mm6*x[3]+1*mm7*x[2]"
,"+8*mm1*((-(-9*x[8+1]*x[8+8]-36*x[8+2]*x[8+7]+72*x[8+1]^2*x[8+7]-84*x[8+3]*x[8+6]+504*x[8+1]*x[8+2]*x[8+6]-504*x[8+1]^3*x[8+6]-126*x[8+4]*x[8+5]+1008*x[8+1]*x[8+3]*x[8+5]+756*
x[8+2]^2*x[8+5]-4536*x[8+1]^2*x[8+2]*x[8+5]+3024*x[8+1]^4*x[8+5]+630*x[8+1]*x[8+4]^2+2520*x[8+2]*x[8+3]*x[8+4]-7560*x[8+1]^2*x[8+3]*x[8+4]-11340*x[8+1]*x[8+2]^2*x[8+4]+30240*x[8+1]^3*x[8+2]*x[8+4]
-15120*x[8+1]^5*x[8+4]+560*x[8+3]^3-15120*x[8+1]*x[8+2]*x[8+3]^2+20160*x[8+1]^3*x[8+3]^2-7560*x[8+2]^3*x[8+3]+90720*x[8+1]^2*x[8+2]^2*x[8+3]-151200*x[8+1]^4*x[8+2]*x[8+3]+60480*
x[8+1]^6*x[8+3]+22680*x[8+1]*x[8+2]^4-151200*x[8+1]^3*x[8+2]^3+272160*x[8+1]^5*x[8+2]^2-181440*x[8+1]^7*x[8+2]+40320*x[8+1]^9))) +28*mm2*x[8]+56*mm3*x[7]+70*mm4*x[6]+56*mm5*x[5]+28*mm6*x[4]+8*mm7*x[3]+1*mm8*x[2]")
if((!missing(beta))&&(length(beta)==2))
{
a..=beta[1]
b..=beta[2]
JA0 =
c("(b..*x[1]+a..)*mm1",
"(b..^2*x[2]+2*a..*b..*x[1]+a..^2)*mm2+(2*b..*x[2]+2*a..*x[1])*mm1",
"(b..^3*x[3]+3*a..*b..^2*x[2]+3*a..^2*b..*x[1]+a..^3)*mm3+(3*b..^2*x[3]+6*a..*b..*x[2]+3*a..^2*x[1])*mm2+(3*b..*x[3]+3*a..*x[2])*mm1",
"(b..^4*x[4]+4*a..*b..^3*x[3]+6*a..^2*b..^2*x[2]+4*a..^3*b..*x[1]+a..^4)*mm4+(4*b..^3*x[4]+12*a..*b..^2*x[3]+12*a..^2*b..*x[2]+4*a..^3*x[1])*mm3+(6*b..^2*x[4]+12*a..*b..*x[3]+6*a..^2*x[2])*mm2+(4*b..*x[4]+4*a..*x[3])*mm1",
"(b..^5*x[5]*mm5+5*a..*b..^4*x[4]*mm5+10*a..^2*b..^3*x[3]*mm5+10*a..^3*b..^2*x[2]*mm5+5*a..^4*b..*x[1]*mm5+a..^5*mm5+5*b..^4*x[5]*mm4+20*a..*b..^3*x[4]*mm4+30*a..^2*b..^2*x[3]*mm4+20*a..^3*b..*x[2]*mm4+5*a..^4*x[1]*mm4+10*b..^3*x[5]*mm3+30*a..*b..^2*x[4]*mm3+30*a..^2*b..*x[3]*mm3+10*a..^3*x[2]*mm3+10*b..^2*x[5]*mm2+20*a..*b..*x[4]*mm2+10*a..^2*x[3]*mm2+5*b..*x[5]*mm1+5*a..*x[4]*mm1)",
"(b..^6*x[6]*mm6+6*a..*b..^5*x[5]*mm6+15*a..^2*b..^4*x[4]*mm6+20*a..^3*b..^3*x[3]*mm6+15*a..^4*b..^2*x[2]*mm6+6*a..^5*b..*x[1]*mm6+a..^6*mm6+6*b..^5*x[6]*mm5+30*a..*b..^4*x[5]*mm5+60*a..^2*b..^3*x[4]*mm5+60*a..^3*b..^2*x[3]*mm5
+30*a..^4*b..*x[2]*mm5+6*a..^5*x[1]*mm5+15*b..^4*x[6]*mm4+60*a..*b..^3*x[5]*mm4+90*a..^2*b..^2*x[4]*mm4+60*a..^3*b..*x[3]*mm4+15*a..^4*x[2]*mm4+20*b..^3*x[6]*mm3+60*a..*b..^2*x[5]*mm3+60*a..^2*b..*x[4]*mm3+20*a..^3*x[3]*mm3+15*
b..^2*x[6]*mm2+30*a..*b..*x[5]*mm2+15*a..^2*x[4]*mm2+6*b..*x[6]*mm1+6*a..*x[5]*mm1)",
"(b..^7*x[7]*mm7+7*a..*b..^6*x[6]*mm7+21*a..^2*b..^5*x[5]*mm7+35*a..^3*b..^4*x[4]*mm7+35*a..^4*b..^3*x[3]*mm7+21*a..^5*b..^2*x[2]*mm7+7*a..^6*b..*x[1]*mm7+a..^7*mm7+7*b..^6*x[7]*mm6+42*a..*b..^5*x[6]*mm6+105*a..^2*b..^4*x[5]*
mm6+140*a..^3*b..^3*x[4]*mm6+105*a..^4*b..^2*x[3]*mm6+42*a..^5*b..*x[2]*mm6+7*a..^6*x[1]*mm6+21*b..^5*x[7]*mm5+105*a..*b..^4*x[6]*mm5+210*a..^2*b..^3*x[5]*mm5+210*a..^3*b..^2*x[4]*mm5+105*a..^4*b..*x[3]*mm5+21*a..^5*x[2]*mm5+
35*b..^4*x[7]*mm4+140*a..*b..^3*x[6]*mm4+210*a..^2*b..^2*x[5]*mm4+140*a..^3*b..*x[4]*mm4+35*a..^4*x[3]*mm4+35*b..^3*x[7]*mm3+105*a..*b..^2*x[6]*mm3+105*a..^2*b..*x[5]*mm3+35*a..^3*x[4]*mm3+21*b..^2*x[7]*mm2+42*a..*b..*x[6]*
mm2+21*a..^2*x[5]*mm2+7*b..*x[7]*mm1+7*a..*x[6]*mm1)",
"(b..^8*x[8]*mm8+8*a..*b..^7*x[7]*mm8+28*a..^2*b..^6*x[6]*mm8+56*a..^3*b..^5*x[5]*mm8+70*a..^4*b..^4*x[4]*mm8+56*a..^5*b..^3*x[3]*mm8+28*a..^6*b..^2*x[2]*mm8+8*a..^7*b..*x[1]*mm8+a..^8*mm8+8*b..^7*x[8]*mm7+56*a..*b..^6*x[7]*mm7
+168*a..^2*b..^5*x[6]*mm7+280*a..^3*b..^4*x[5]*mm7+280*a..^4*b..^3*x[4]*mm7+168*a..^5*b..^2*x[3]*mm7+56*a..^6*b..*x[2]*mm7+8*a..^7*x[1]*mm7+28*b..^6*x[8]*mm6+168*a..*b..^5*x[7]*mm6+420*a..^2*b..^4*x[6]*mm6+560*a..^3*b..^3*x[5]*mm6
+420*a..^4*b..^2*x[4]*mm6+168*a..^5*b..*x[3]*mm6+28*a..^6*x[2]*mm6+56*b..^5*x[8]*mm5+280*a..*b..^4*x[7]*mm5+560*a..^2*b..^3*x[6]*mm5+560*a..^3*b..^2*x[5]*mm5+280*a..^4*b..*x[4]*mm5+56*a..^5*x[3]*mm5+70*b..^4*x[8]*mm4+280*
a..*b..^3*x[7]*mm4+420*a..^2*b..^2*x[6]*mm4+280*a..^3*b..*x[5]*mm4+70*a..^4*x[4]*mm4+56*b..^3*x[8]*mm3+168*a..*b..^2*x[7]*mm3+168*a..^2*b..*x[6]*mm3+56*a..^3*x[5]*mm3+28*b..^2*x[8]*mm2+56*a..*b..*x[7]*mm2+28*a..^2*x[6]*mm2+8
*b..*x[8]*mm1+8*a..*x[7]*mm1)")
print('Invoked special jump structure.')
}
if(checklist[7]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[7])[2],')*(',JA0[i],')')
}
}
if(checklist[8]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[8])[2],')*(',JA1[i],')')
}
}
if(checklist[9]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[9])[2],')*(',JA2[i],')')
}
}
}
#c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta')
if(checklist[7]==0){Lam0=function(t){0}}
if(checklist[8]==0){Lam1=function(t){0}}
if(checklist[9]==0){Lam2=function(t){0}}
if(checklist[10]==0){Jmu=function(t){0}}
if(checklist[11]==0){Jsig=function(t){0}}
if(checklist[12]==0){Jlam=function(t){0}}
if(checklist[13]==0){Jalpha=function(t){0}}
if(checklist[14]==0){Jbeta=function(t){0}}
if(checklist[15]==0){Ja=function(t){0}}
if(checklist[16]==0){Jb=function(t){0}}
odekernel=paste0('c(',dims[1],'\n,',dims[2],'\n,',dims[3],'\n,',dims[4],'\n,',dims[5],'\n,',dims[6],'\n,',dims[7],'\n,',dims[8],'\n,',
dims2[1],'\n,',dims2[2],'\n,',dims2[3],'\n,',dims2[4],'\n,',dims2[5],'\n,',dims2[6],'\n,',dims2[7],'\n,',dims2[8],'\n,',
body(namess[7])[2],'\n,',paste('(',body(namess[8])[2],')*x[9]'),'\n,',paste('(',body(namess[9])[2],')*x[10]'),')')
}
if(TR.order ==4)
{
MAT=rbind(
c('1','1*x[1]','1*x[2]','','',''),
c('2*x[1]','2*x[2]','2*x[3]','1','1*x[1]','1*x[2]'),
c('3*x[2]','3*x[3]','3*x[4]','3*x[1]','3*x[2]','3*x[3]'),
c('4*x[3]','4*x[4]','4*x[5]','6*x[2]','6*x[3]','6*x[4]'))
MAT2=rbind(
c('1','1*x[4+1]','1*x[4+2]','','',''),
c('2*x[4+1]','2*x[4+2]','2*x[4+3]','1','1*x[4+1]','1*x[4+2]'),
c('3*x[4+2]','3*x[4+3]','3*x[4+4]','3*x[4+1]','3*x[4+2]','3*x[4+3]'),
c('4*x[4+3]','4*x[4+4]','4*x[4+5]','6*x[4+2]','6*x[4+3]','6*x[4+4]'))
dims=rep('',4)
dims2=rep('',4)
for(i in which(checklist[1:6]==1))
{
for(j in 1:4)
{
if(MAT[j,i]!='')
{
dims[j] = paste0(dims[j],'+(',body(namess[i])[2],')*',MAT[j,i])
dims2[j] = paste0(dims2[j],'+(',body(namess[i])[2],')*',MAT2[j,i])
}
}
}
if(Jtype=='Add')
{
JA0 =
c("+1*mm1"
,"+2*mm1*x[1]+1*mm2"
,"+3*mm1*x[2]+3*mm2*x[1]+1*mm3"
,"+4*mm1*x[3]+6*mm2*x[2]+4*mm3*x[1]+1*mm4")
JA1 =
c("+1*mm1*x[1]"
,"+2*mm1*x[2]+1*mm2*x[1]"
,"+3*mm1*x[3]+3*mm2*x[2]+1*mm3*x[1]"
,"+4*mm1*x[4]+6*mm2*x[3]+4*mm3*x[2]+1*mm4*x[1]")
JA2 =
c("+1*mm1*x[2]"
,"+2*mm1*x[3]+1*mm2*x[2]"
,"+3*mm1*x[4]+3*mm2*x[3]+1*mm3*x[2]"
,"+4*mm1*x[5]+6*mm2*x[4]+4*mm3*x[3]+1*mm4*x[2]")
}
if((!missing(beta))&&(length(beta)==2))
{
a..=beta[1]
b..=beta[2]
JA0 =
c("b..*x[1]*mm1+a..*mm1"
,"b..^2*x[2]*mm2+2*a..*b..*x[1]*mm2+a..^2*mm2+2*b..*x[2]*mm1+2*a..*x[1]*mm1"
,"b..^3*x[3]*mm3+3*a..*b..^2*x[2]*mm3+3*a..^2*b..*x[1]*mm3+a..^3*mm3+3*b..^2*x[3]*mm2+6*a..*b..*x[2]*mm2+3*a..^2*x[1]*mm2+3*b..*x[3]*mm1+3*a..*x[2]*mm1"
,"b..^4*x[4]*mm4+4*a..*b..^3*x[3]*mm4+6*a..^2*b..^2*x[2]*mm4+4*a..^3*b..*x[1]*mm4+a..^4*mm4+4*b..^3*x[4]*mm3+12*a..*b..^2*x[3]*mm3+12*a..^2*b..*x[2]*mm3+4*a..^3*x[1]*mm3+6*b..^2*x[4]*mm2+12*a..*b..*x[3]*mm2+6*a..^2*x[2]*mm2+4*b..*x[4]*mm1+4*a..*x[3]*mm1")
print('Invoked special jump structure.')
}
if(checklist[7]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[7])[2],')*(',JA0[i],')')
}
}
if(checklist[8]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[8])[2],')*(',JA1[i],')')
}
}
if(checklist[9]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[9])[2],')*(',JA2[i],')')
}
}
if(checklist[7]==0)
{
Lam0=function(t){0}
}
if(checklist[8]==0)
{
Lam1=function(t){0}
}
if(checklist[9]==0)
{
Lam2=function(t){0}
}
odekernel=paste0('c(',dims[1],'\n,',dims[2],'\n,',dims[3],'\n,',dims[4],'\n,',
dims2[1],'\n,',dims2[2],'\n,',dims2[3],'\n,',dims2[4],'\n,',
body(namess[7])[2],'\n,',paste('(',body(namess[8])[2],')*x[5]'),'\n,',paste('(',body(namess[9])[2],')*x[6]'),')')
}
if(Jdist=='Normal')
{
prem =
'mm1 = mu
mm2 = (mu^2)+ (sig^2)
mm3 = (mu^3)+ 3* (mu^1)*(sig^2)
mm4 = (mu^4)+ 6* (mu^2)*(sig^2)+3*(sig^4)
mm5 = (mu^5)+ 10*(mu^3)*(sig^2)+15*(mu^1)*(sig^4)
mm6 = (mu^6)+ 15*(mu^4)*(sig^2)+45*(mu^2)*(sig^4)+15*(sig^6)
mm7 = (mu^7)+ 21*(mu^5)*(sig^2)+105*(mu^3)*(sig^4)+105*mu*(sig^6)
mm8 = (mu^8)+ 28*(mu^6)*(sig^2)+210*(mu^4)*(sig^4)+420*(mu^2)*(sig^6)+105*(sig^8)
'
prem=paste0('mu = ',body('Jmu')[2],'\n','sig = ', body('Jsig')[2],'\n',prem)
}
if(Jdist=='Exponential')
{
prem =
'mm1=1*(mua^1)
mm2=2*(mua^2)
mm3=6*(mua^3)
mm4=24*(mua^4)
mm5=120*(mua^5)
mm6=720*(mua^6)
mm7=5040*(mua^7)
mm8=40320*(mua^8)
'
prem=paste0('mua = ',body('Jlam')[2],'\n',prem)
}
if(Jdist=='Gamma')
{
prem =
'mm1=(alph*bet)
mm2=(alph+1)*alph*bet^2
mm3=(alph+2)*(alph+1)*alph*bet^3
mm4=(alph+3)*(alph+2)*(alph+1)*alph*bet^4
mm5=(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^5
mm6=(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^6
mm7=(alph+6)*(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^7
mm8=(alph+7)*(alph+6)*(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^8
'
prem=paste0('alph = ',body('Jalpha')[2],'\n',prem)
prem=paste0('bet = ',body('Jbeta')[2],'\n',prem)
}
if(Jdist=='Laplace')
{
#print('yeah')
prem =
'mm1=0.5*(+2*aaaa^1*bbbb^0)
mm2=0.5*(+2*aaaa^2*bbbb^0+4*aaaa^0*bbbb^2)
mm3=0.5*(+2*aaaa^3*bbbb^0+12*aaaa^1*bbbb^2)
mm4=0.5*(+2*aaaa^4*bbbb^0+24*aaaa^2*bbbb^2+48*aaaa^0*bbbb^4)
mm5=0.5*(+2*aaaa^5*bbbb^0+40*aaaa^3*bbbb^2+240*aaaa^1*bbbb^4)
mm6=0.5*(+2*aaaa^6*bbbb^0+60*aaaa^4*bbbb^2+720*aaaa^2*bbbb^4+1440*aaaa^0*bbbb^6)
mm7=0.5*(+2*aaaa^7*bbbb^0+84*aaaa^5*bbbb^2+1680*aaaa^3*bbbb^4+10080*aaaa^1*bbbb^6)
mm8=0.5*(+2*aaaa^8*bbbb^0+112*aaaa^6*bbbb^2+3360*aaaa^4*bbbb^4+40320*aaaa^2*bbbb^6+80640*aaaa^0*bbbb^8)
'
prem=paste0('aaaa = ',body('Ja')[2],'\n',prem)
prem=paste0('bbbb = ',body('Jb')[2],'\n',prem)
}
if(Jdist=='Uniform')
{
prem =
'
mm1=(bbbb+aaaa)/2
mm2=(aaaa^2+aaaa*bbbb +bbbb^2)/3
mm3=1/4*(aaaa+bbbb)*(aaaa^2+bbbb^2)
mm4=1/5*(aaaa^4+aaaa^3*bbbb+aaaa^2*bbbb^2+aaaa*bbbb^3+bbbb^4)
mm5=1/6*(aaaa^5+aaaa^4*bbbb+aaaa^3*bbbb^2+aaaa^2*bbbb^3+aaaa^1*bbbb^4+bbbb^5)
mm6=1/7*(aaaa^6+aaaa^5*bbbb^1+aaaa^4*bbbb^2+aaaa^3*bbbb^3+aaaa^2*bbbb^4+aaaa^1*bbbb^5+bbbb^6)
mm7=1/8*(aaaa^7+aaaa^6*bbbb^1+aaaa^5*bbbb^2+aaaa^4*bbbb^3+aaaa^3*bbbb^4+aaaa^2*bbbb^5+aaaa^1*bbbb^6+bbbb^7)
mm8=1/9*(aaaa^8+aaaa^7*bbbb^1+aaaa^6*bbbb^2+aaaa^5*bbbb^3+aaaa^4*bbbb^4+aaaa^3*bbbb^5+aaaa^2*bbbb^6+aaaa^1*bbbb^7+bbbb^8)
'
prem=paste0('aaaa = ',body('Ja')[2],'\n',prem)
prem=paste0('bbbb = ',body('Jb')[2],'\n',prem)
}
odekernel=paste0('{',prem,odekernel,'}')
write(odekernel,'res.txt')
ff <- function(x,t){}
body(ff) = (parse(text=odekernel))
t.alpha=
c(0.000000000000000000000000000000000000000000000000000000000000,
0.100000000000000000000000000000000000000000000000000000000000,
0.539357840802981787532485197881302436857273449701009015505500,
0.809036761204472681298727796821953655285910174551513523258250,
0.309036761204472681298727796821953655285910174551513523258250,
0.981074190219795268254879548310562080489056746118724882027805,
0.833333333333333333333333333333333333333333333333333333333333,
0.354017365856802376329264185948796742115824053807373968324184,
0.882527661964732346425501486979669075182867844268052119663791,
0.642615758240322548157075497020439535959501736363212695909875,
0.357384241759677451842924502979560464040498263636787304090125,
0.117472338035267653574498513020330924817132155731947880336209,
0.833333333333333333333333333333333333333333333333333333333333,
0.309036761204472681298727796821953655285910174551513523258250,
0.539357840802981787532485197881302436857273449701009015505500,
0.100000000000000000000000000000000000000000000000000000000000,
1.00000000000000000000000000000000000000000000000000000000000)
N=round((t-s)/delt)
ttt=seq(s,t,by=delt)
solver=function()
{
MM=matrix(0,2*TR.order+3,N)
MA= c(rep(Xs^(1:TR.order),2),0,0,0)
tt=s
MM[,1]=MA
ii=0
while(ii < N)
{
ii=ii+1
x0=MA
fx0=ff(x0,tt)
x1=x0+delt*(0.1*fx0)
fx1=ff(x1,tt+t.alpha[2]*delt)
x2=x0+delt*(-0.915176561375291*fx0 +1.45453440217827*fx1)
fx2=ff(x2,tt+t.alpha[3]*delt)
x3=x0+delt*( 0.202259190301118*fx0 +0.606777570903354*fx2)
fx3=ff(x3,tt+t.alpha[4]*delt)
x4=x0+delt*( 0.184024714708644*fx0 +0.197966831227192*fx2-0.0729547847313633*fx3)
fx4=ff(x4,tt+t.alpha[5]*delt)
x5=x0+delt*( 0.0879007340206681*fx0 +0.410459702520261*fx3+0.482713753678866*fx4)
fx5=ff(x5,tt+t.alpha[6]*delt)
x6=x0+delt*(0.085970050490246*fx0 +0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5)
fx6=ff(x6,tt+t.alpha[7]*delt)
x7=x0+delt*(0.120930449125334*fx0 +0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6)
fx7=ff(x7,tt+t.alpha[8]*delt)
x8=x0+delt*(0.110854379580391*fx0 -0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7)
fx8=ff(x8,tt+t.alpha[9]*delt)
x9=x0+delt*(0.112054414752879*fx0 -0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8)
fx9=ff(x9,tt+t.alpha[10]*delt)
x10=x0+delt*(0.113976783964186*fx0 -0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9)
fx10=ff(x10,tt+t.alpha[11]*delt)
x11=x0+delt*(0.0798314528280196*fx0 -0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10)
fx11=ff(x11,tt+t.alpha[12]*delt)
x12=x0+delt*(0.985115610164857*fx0 +0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11)
fx12=ff(x12,tt+t.alpha[13]*delt)
x13=x0+delt*(0.895080295771633*fx0 +0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12)
fx13=ff(x13,tt+t.alpha[14]*delt)
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13)
fx14=ff(x14,tt+t.alpha[15]*delt)
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14)
fx15=ff(x15,tt+t.alpha[16]*delt)
x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15)
fx16=ff(x16,tt+t.alpha[17]*delt)
MA=MA+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt
tt=tt+delt
MM[,ii]=MA
setTxtProgressBar(pb,ii," "," ")
}
return(MM)
}
pb <- txtProgressBar(1,N,1,style = 1,width = 65)
MM=solver()
close(pb)
if(eval.density)
{
DD=matrix(0,length(Xt),N)
if(Dtype=='Saddle.Gamma')
{
print('yeaahh')
dens1=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
dens2=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p0 = exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
p1 = 1-p0
VVV[1,] = (MM[1,]-p0*MM[1+TR.order,])/p1
VVV[2,] = (MM[2,]-p0*MM[2+TR.order,])/p1
VVV[3,] = (MM[3,]-p0*MM[3+TR.order,])/p1
VVV[4,] = (MM[4,]-p0*MM[4+TR.order,])/p1
if(!factorize)
{
for(i in 1:N)
{
DD[,i]=dens1(Xt,MM[1:4,i])
}
}
if(factorize)
{
if(factor.type=='Diffusion')
{
for(i in 1:N)
{
if(p0[i]>0.01){DD[,i]=dens1(Xt,MMM[1:4,i])*p0[i]+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
if(factor.type=='Hawkes')
{
if(p0[i]>0.01){DD[,i]=p0[i]/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
DD=list(density=DD,MSH=0)
}
}
if(Dtype=='Gamma.Saddle')
{
print('yeahh')
dens2=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
dens1=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p0 = exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
p1 = 1-p0
VVV[1,] = (MM[1,]-p0*MM[1+TR.order,])/p1
VVV[2,] = (MM[2,]-p0*MM[2+TR.order,])/p1
VVV[3,] = (MM[3,]-p0*MM[3+TR.order,])/p1
VVV[4,] = (MM[4,]-p0*MM[4+TR.order,])/p1
if(!factorize)
{
for(i in 1:N)
{
DD[,i]=dens1(Xt,MM[1:4,i])
}
}
if(factorize)
{
if(factor.type=='Diffusion')
{
for(i in 1:N)
{
if(p0[i]>0.01){DD[,i]=dens1(Xt,MMM[1:4,i])*p0[i]+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
if(factor.type=='Hawkes')
{
if(p0[i]>0.01){DD[,i]=p0[i]/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
DD=list(density=DD,MSH=0)
}
}
if(Dtype=='Gamma.A')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=dnorm(Xt,MMM[1,i],sqrt(MMM[2,i]-MMM[1,i]^2))*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Normal.A')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=dnorm(Xt,MMM[1,i],sqrt(MMM[2,i]-MMM[1,i]^2))*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
return(dnorm(xx,k[1],sqrt(k[2])))
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Saddlepoint')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=saddlep(Xt,MMM[1:4,i])*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])*Xt*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Edgeworth')
{
edgeworth=function(x,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
v3=k[3]/(sqrt(k[2])^3)
v4=k[4]/(sqrt(k[2])^4)
z=(x-k[1])/sqrt(k[2])
h2=function(z){z^2-1}
h3=function(z){z^3-3*z}
h4=function(z){z^4-6*z^2+3}
h5=function(z){z^5-10*z^3+15*z}
return(pmax(dnorm(x,k[1],sqrt(k[2]))*(1+v3/6*h3(z)+v4/24*h4(z)),0))
}
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=edgeworth(Xt,MMM[1:4,i])*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+edgeworth(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
VVV[1,] = MM[1,]#-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-p*MM[1+TR.order,]^2)/(1-p)
VVV[3,] = (MM[3,]-p*MM[1+TR.order,]^3)/(1-p)
VVV[4,] = (MM[4,]-p*MM[1+TR.order,]^4)/(1-p)
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])*Xt*(abs(Xt-MM[1,i])<0.005)+edgeworth(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
for(i in 1:N)
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order Edgeworth not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order Edgeworth not supported at present!')
}
}
if(DTR.order==4)
{
pearson=function(x,k,type='Normal')
{
dx=diff(x)[1]
u = k*0
u[1,]=k[1,]
u[2,]=k[2,]
u[3,]=k[3,]
u[4,]=k[4,]
det=1*u[2,]*u[4,]+u[1,]*u[3,]*u[2,]+u[2,]*u[1,]*u[3,]-1*u[3,]*u[3,]-u[2,]*u[2,]*u[2,]-u[1,]*u[1,]*u[4,]
b11=(u[2,]*u[4,]-u[3,]*u[3,])/det
b12=(u[2,]*u[3,]-u[1,]*u[4,])/det
b13=(u[1,]*u[3,]-u[2,]*u[2,])/det
b21=(u[3,]*u[2,]-u[1,]*u[4,])/det
b22=(1*u[4,]-u[2,]*u[2,])/det
b23=(u[2,]*u[1,]-1*u[3,])/det
b31=(u[1,]*u[3,]-u[2,]*u[2,])/det
b32=(u[1,]*u[2,]-1*u[3,])/det
b33=(1*u[2,]-u[1,]*u[1,])/det
#Ntype
if(type=='Normal')
{
V1=0
V2=1
V3=2*u[1,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
thetas1=-betas1/1
thetas2=-betas2/2
thetas3=-betas3/3
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum(exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3)*rho*diff(y)[1])
DD[,i] = exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3)/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Gtype
if(type=='Gamma')
{
V1=1
V2=2*u[1,]
V3=3*u[2,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alphas =1-betas1
thetas1=-betas2/1
thetas2=-betas3/2
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum(tau^(alphas[i]-1)*exp(thetas1[i]*tau+thetas2[i]*tau^2)*rho*diff(y)[1])
DD[,i] = (x^(alphas[i]-1)*exp(thetas1[i]*x+thetas2[i]*x^2))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# IGtype
if(type=='InvGamma')
{
V1=2*u[1,]
V2=3*u[2,]
V3=4*u[3,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alphas =-betas2
nu =betas1
thetas =-betas3
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum((tau^alphas[i])*exp(nu[i]/tau+thetas[i]*tau)*rho*diff(y)[1])
DD[,i] = (x^alphas[i])*exp(nu[i]/x+thetas[i]*x)/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Btype
if(type=='Beta')
{
V1=1-2*u[1,]
V2=2*u[1,]-3*u[2,]
V3=3*u[2,]-4*u[3,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alpha =1-betas1
nu =1+(betas1+betas2+betas3)
theta1 =betas3
DD=matrix(0,length(x),length(u[1,]))
for(i in 1:length(x))
{
DD[i,]=(x[i]^(alpha-1))*((1-x[i])^(nu-1))*exp(theta1*x[i])
}
for(i in 1:dim(u)[2])
{
DD[,i]=DD[,i]/sum(DD[,i]*dx)
}
return(list(density=DD,MSH=0))
}
}
}
if(DTR.order==8)
{
pearson=function(x,k,type='Normal')
{
dx=diff(x)[1]
u1= k[1,]
u2=k[2,]
u3=k[3,]
u4=k[4,]
u5=k[5,]
u6=k[6,]
u7=k[7,]
u8=k[8,]
det=-u1*(u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*
(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))-u5*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5)))+u2*(u1*
(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*
(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4)))-u3*(u1*
(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))-u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u3*
(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))+u2*
(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))+u4*
(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u4*(u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*
(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*
(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))-u5*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))
b11= (u2*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))+u4*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))-u5*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5)))/det
b12= (-u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))+u2*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u3*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u4*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5)))/det
b13= (u1*(u3*(u6*u8-u7*u7)-u4*(u5*u8-u6*u7)+u5*(u5*u7-u6*u6))-u2*(u2*(u6*u8-u7*u7)-u4*(u4*u8-u5*u7)+u5*(u4*u7-u5*u6))+u3*(u2*(u5*u8-u6*u7)-u3*(u4*u8-u5*u7)+u5*(u4*u6-u5*u5))-u4*(u2*(u5*u7-u6*u6)-u3*(u4*u7-u5*u6)+u4*(u4*u6-u5*u5)))/det
b14= (-u1*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u6*u6)+u5*(u4*u7-u5*u6))+u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u5*u6)+u5*(u3*u7-u5*u5))-u3*(u2*(u4*u8-u6*u6)-u3*(u3*u8-u5*u6)+u5*(u3*u6-u4*u5))+u4*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u5*u5)+u4*(u3*u6-u4*u5)))/det
b15= (u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b21= (-u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))+u3*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))-u4*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u5*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5)))/det
b22= (-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u3*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)-u4*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u6*(u5*u7-u6*u6))/det
b23= (u2*(u1*(u6*u8-u7*u7)-u4*(u3*u8-u4*u7)+u5*(u3*u7-u4*u6))-u3*(u1*(u5*u8-u6*u7)-u3*(u3*u8-u4*u7)+u5*(u3*u6-u4*u5))-u3*(u6*u8-u7*u7)+u4*(u5*u8-u6*u7)+u4*(u1*(u5*u7-u6*u6)-u3*(u3*u7-u4*u6)+u4*(u3*u6-u4*u5))-u5*(u5*u7-u6*u6))/det
b24= (-u2*(u1*(u5*u8-u6*u7)-u4*(u2*u8-u4*u6)+u5*(u2*u7-u4*u5))+u3*(u1*(u4*u8-u6*u6)-u3*(u2*u8-u4*u6)+u5*(u2*u6-u4*u4))+u3*(u5*u8-u6*u7)-u4*(u4*u8-u6*u6)-u4*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u4*u5)+u4*(u2*u6-u4*u4))+u5*(u4*u7-u5*u6))/det
b25= (u2*(u1*(u5*u7-u6*u6)-u4*(u2*u7-u3*u6)+u5*(u2*u6-u3*u5))-u3*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u3*u6)+u5*(u2*u5-u3*u4))-u3*(u5*u7-u6*u6)+u4*(u4*u7-u5*u6)+u4*(u1*(u4*u6-u5*u5)-u3*(u2*u6-u3*u5)+u4*(u2*u5-u3*u4))-u5*(u4*u6-u5*u5))/det
b31= (u1*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4)))/det
b32= (u1*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))-u3*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u3*(u6*u8-u7*u7)+u5*(u4*u8-u5*u7)+u4*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u6*(u4*u7-u5*u6))/det
b33= (-u1*(u1*(u6*u8-u7*u7)-u4*(u3*u8-u4*u7)+u5*(u3*u7-u4*u6))+u3*(u1*(u4*u8-u5*u7)-u2*(u3*u8-u4*u7)+u5*(u3*u5-u4*u4))+u2*(u6*u8-u7*u7)-u4*(u4*u8-u5*u7)-u4*(u1*(u4*u7-u5*u6)-u2*(u3*u7-u4*u6)+u4*(u3*u5-u4*u4))+u5*(u4*u7-u5*u6))/det
b34= (u1*(u1*(u5*u8-u6*u7)-u4*(u2*u8-u4*u6)+u5*(u2*u7-u4*u5))-u3*(u1*(u3*u8-u5*u6)-u2*(u2*u8-u4*u6)+u5*(u2*u5-u3*u4))-u2*(u5*u8-u6*u7)+u4*(u3*u8-u5*u6)+u4*(u1*(u3*u7-u5*u5)-u2*(u2*u7-u4*u5)+u4*(u2*u5-u3*u4))-u5*(u3*u7-u5*u5))/det
b35= (-u1*(u1*(u5*u7-u6*u6)-u4*(u2*u7-u3*u6)+u5*(u2*u6-u3*u5))+u3*(u1*(u3*u7-u4*u6)-u2*(u2*u7-u3*u6)+(u2*u4-u3*u3)*u5)+u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)-u4*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u3*u5)+u4*(u2*u4-u3*u3))+u5*(u3*u6-u4*u5))/det
b41= (-u1*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))-u3*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)+u5*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b42= (-u1*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u2*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)+u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4))+u6*(u4*u6-u5*u5))/det
b43= (u1*(u1*(u5*u8-u6*u7)-u3*(u3*u8-u4*u7)+u5*(u3*u6-u4*u5))-u2*(u1*(u4*u8-u5*u7)-u2*(u3*u8-u4*u7)+u5*(u3*u5-u4*u4))-u2*(u5*u8-u6*u7)+u3*(u4*u8-u5*u7)+u4*(u1*(u4*u6-u5*u5)-u2*(u3*u6-u4*u5)+u3*(u3*u5-u4*u4))-u5*(u4*u6-u5*u5))/det
b44= (-u1*(u1*(u4*u8-u6*u6)-u3*(u2*u8-u4*u6)+u5*(u2*u6-u4*u4))+u2*(u1*(u3*u8-u5*u6)-u2*(u2*u8-u4*u6)+u5*(u2*u5-u3*u4))+u2*(u4*u8-u6*u6)-u3*(u3*u8-u5*u6)-u4*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u4*u4)+u3*(u2*u5-u3*u4))+u5*(u3*u6-u4*u5))/det
b45= (u1*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u3*u6)+u5*(u2*u5-u3*u4))-u2*(u1*(u3*u7-u4*u6)-u2*(u2*u7-u3*u6)+(u2*u4-u3*u3)*u5)-u2*(u4*u7-u5*u6)+u3*(u3*u7-u4*u6)+u4*(u1*(u3*u5-u4*u4)-u2*(u2*u5-u3*u4)+u3*(u2*u4-u3*u3))-u5*(u3*u5-u4*u4))/det
b51= (u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b52= (u1*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))-u2*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u3*(u5*u7-u6*u6)+u4*(u4*u7-u5*u6)+u3*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4))-u5*(u4*u6-u5*u5))/det
b53= (-u1*(u1*(u5*u7-u6*u6)-u3*(u3*u7-u4*u6)+u4*(u3*u6-u4*u5))+u2*(u1*(u4*u7-u5*u6)-u2*(u3*u7-u4*u6)+u4*(u3*u5-u4*u4))+u2*(u5*u7-u6*u6)-u3*(u4*u7-u5*u6)-u3*(u1*(u4*u6-u5*u5)-u2*(u3*u6-u4*u5)+u3*(u3*u5-u4*u4))+u4*(u4*u6-u5*u5))/det
b54= (u1*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u4*u5)+u4*(u2*u6-u4*u4))-u2*(u1*(u3*u7-u5*u5)-u2*(u2*u7-u4*u5)+u4*(u2*u5-u3*u4))-u2*(u4*u7-u5*u6)+u3*(u3*u7-u5*u5)+u3*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u4*u4)+u3*(u2*u5-u3*u4))-u4*(u3*u6-u4*u5))/det
b55= (-u1*(u1*(u4*u6-u5*u5)-u3*(u2*u6-u3*u5)+u4*(u2*u5-u3*u4))+u2*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u3*u5)+u4*(u2*u4-u3*u3))+u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)-u3*(u1*(u3*u5-u4*u4)-u2*(u2*u5-u3*u4)+u3*(u2*u4-u3*u3))+u4*(u3*u5-u4*u4))/det
#Ntype
if(type=='Normal')
{
V1=0
V2=1
V3=2*u1
V4=3*u2
V5=4*u3
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
thetas1=-betas1/1
thetas2=-betas2/2
thetas3=-betas3/3
thetas4=-betas4/4
thetas5=-betas5/5
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum(exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3+thetas4[i]*tau^4+thetas5[i]*tau^5)*rho*diff(y)[1])
DD[,i] = (exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3+thetas4[i]*x^4+thetas5[i]*x^5))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Gtype
if(type=='Gamma')
{
V1=1
V2=2*u1
V3=3*u2
V4=4*u3
V5=5*u4
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alphas =1-betas1
thetas1=-betas2/1
thetas2=-betas3/2
thetas3=-betas4/3
thetas4=-betas5/4
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum(tau^(alphas[i]-1)*exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3+thetas4[i]*tau^4)*rho*diff(y)[1])
DD[,i] = (x^(alphas[i]-1)*exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3+thetas4[i]*x^4))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# IGtype
if(type=='InvGamma')
{
V1=2*u1
V2=3*u2
V3=4*u3
V4=5*u4
V5=6*u5
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alphas =-betas2
nu =betas1
theta1 =-betas3
theta2 =-betas4/2
theta3 =-betas5/3
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum((tau^alphas[i])*exp(nu[i]/tau+theta1[i]*tau+theta2[i]*tau^2+theta3[i]*tau^3)*rho*diff(y)[1])
DD[,i] = ((x^alphas[i])*exp(nu[i]/x+theta1[i]*x+theta2[i]*x^2+theta3[i]*x^3))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Btype
if(type=='Beta')
{
V1=1-2*u1
V2=2*u1-3*u2
V3=3*u2-4*u3
V4=4*u3-5*u4
V5=5*u4-6*u5
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alpha =1-betas1
nu =1+(betas1+betas2+betas3+betas4+betas5)
theta1 =betas3+betas4+betas5
theta2 =(betas4+betas5)/2
theta3 =betas5/3
DD=matrix(0,length(x),length(u1))
for(i in 1:length(x))
{
DD[i,]=(x[i]^(alpha-1))*((1-x[i])^(nu-1))*exp(theta1*x[i]+theta2*x[i]^2+theta3*x[i]^3)
}
for(i in 1:length(u1))
{
DD[,i]=DD[,i]/sum(DD[,i]*dx)
}
return(list(density=DD,MSH=0))
}
}
}
if(Dtype=='Normal')
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Normal')
#print(d1)
d2=pearson(Xt,VVV,'Normal')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MMM,'Normal')
}
}
if(Dtype=='Gamma')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must be strictly positive for Dtype = Gamma. ')}
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Gamma')
#print(d1)
d2=pearson(Xt,VVV,'Gamma')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MMM,'Gamma')
}
}
if(Dtype=='InvGamma')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must be strictly positive for Dtype = InvGamma. ')}
DD=pearson(Xt,MM,'InvGamma')
}
if(Dtype=='Beta')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must in [0,1] Dtype = Beta. Try seq(delta,1-delta,delta) for delta>0.')}
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(TR.order==8)
{
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
}
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Beta')
#print(d1)
d2=pearson(Xt,VVV,'Beta')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MM,'Beta')
}
}
}
}
if(!eval.density)
{
DD=list(density=NULL,MSH=NULL);
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
}
# Fix this!!!
if(TR.order==8)
{
u=MM[1:8,]
k=u*0
k[1,]= u[1,]
k[2,]= u[2,]-1*k[1,]*u[1,]
k[3,]= u[3,]-1*k[1,]*u[2,]-2*k[2,]*u[1,]
k[4,]= u[4,]-1*k[1,]*u[3,]-3*k[2,]*u[2,]-3*k[3,]*u[1,]
k[5,]= u[5,]-1*k[1,]*u[4,]-4*k[2,]*u[3,]-6*k[3,]*u[2,]-4*k[4,]*u[1,]
k[6,]= u[6,]-1*k[1,]*u[5,]-5*k[2,]*u[4,]-10*k[3,]*u[3,]-10*k[4,]*u[2,]-5*k[5,]*u[1,]
k[7,]= u[7,]-1*k[1,]*u[6,]-6*k[2,]*u[5,]-15*k[3,]*u[4,]-20*k[4,]*u[3,]-15*k[5,]*u[2,]-6*k[6,]*u[1,]
k[8,]= u[8,]-1*k[1,]*u[7,]-7*k[2,]*u[6,]-21*k[3,]*u[5,]-35*k[4,]*u[4,]-35*k[5,]*u[3,]-21*k[6,]*u[2,]-7*k[7,]*u[1,]
}
if(TR.order==4)
{
u=MM[1:4,]
k=u*0
k[1,]= u[1,]
k[2,]= u[2,]-1*k[1,]*u[1,]
k[3,]= u[3,]-1*k[1,]*u[2,]-2*k[2,]*u[1,]
k[4,]= u[4,]-1*k[1,]*u[3,]-3*k[2,]*u[2,]-3*k[3,]*u[1,]
}
ret= list(density=DD$density,Xt=Xt,time=seq(s,t,delt)[-1],cumulants=k,moments=MM,mesh=DD$MSH,zero_jump_prob =exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]),excess_moments =VVV)
class(ret) = 'JGQD.density'
return(ret)
}
eta21/DiffusionRjgqd documentation built on May 16, 2019, 8:54 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
globalVariables(c('mua','alph'))
JGQD.density=function(Xs=4,Xt=seq(5,8,1/10),s=0,t=5,delt=1/100,Jdist='Normal',Jtype='Add',Dtype='Saddlepoint',Trunc=c(8,4),factorize=FALSE,factor.type='Diffusion',beta,print.output=TRUE,eval.density=TRUE)
{
Integ = TRUE
P=100;alpha=0.1;lower=0;upper=50
TR.order=Trunc[1]
DTR.order=Trunc[2]
Dtypes =c('Saddlepoint','Edgeworth','Normal.A','Gamma.A','Beta','Saddle.Gamma','Gamma.Saddle')
check_for_model=function()
{
txt=''
namess=c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta','Ja','Jb')
func.list=rep(0,length(namess))
obs=objects(pos=1)
for(i in 1:length(namess))
{
if(sum(obs==namess[i])){func.list[i]=1}
}
check=F
if(sum(func.list)==0)
{
txt='
--------------------------------------------------------------------------------
No model has been defined yet! Try for example:
--------------------------------------------------------------------------------
GQD.remove()
G0=function(t){5+2*sin(4*pi*t)}
G1=function(t){-1}
Q1=function(t){1}
Lam0 = function(t){0.1}
Jmu = function(t){0}
Jsig = function(t){1}
model=JGQD.density(5,seq(0,10,1/10),0,5,1/100)
JGQD.plot(model)
--------------------------------------------------------------------------------
'
check=T
}
if((sum(func.list)>0)&&(sum(func.list[c(1:3)+3])==0))
{
txt2='
--------------------------------------------------------------------------------
No diffusion coefficients specified. The resulting process will be that of a
Hawkes process.
{factor.type} changed to "Hawkes"
--------------------------------------------------------------------------------
'
warning(txt2)
factor.type = 'Hawkes'
}
return(list(check=check,txt=txt))
}
check_for=check_for_model()
if(check_for[[1]]){stop(check_for[[2]])}
TR.order=Trunc[1]
DTR.order=Trunc[2]
#Dtypes =c('Saddlepoint','Normal','Beta')
Dindex = which(Dtypes==Dtype)
JDtypes =c('Normal','Exponential','Gamma','Laplace','Uniform')
JDindex = which(JDtypes==Jdist)
b1 = '\n==============================================================================\n'
b2 = '==============================================================================\n'
warn=c(
'Input 1: Argument {lower} must be of length 1.\n'
,'Input 2: Argument {upper} must be of length 1.\n'
,'Input 3: {upper} must be > {lower}.\n'
,'Input 4: Range [lower,upper] must be strictly positive for Dtype Gamma or InvGamma.\n'
,'Input 5: Dtype cannot be Beta for observations not in (0,1).\n'
,'Input 6: P must be >= 10.\n'
,'Input 7: Dtype has to be one of Saddlepoint or Edgeworth.\n'
,'Input 8: Trunc[2] must be <= Trunc[1].\n'
,'Input 9: Large {delt} may result in poor approximations.\n'
,'Input 10: Argument {delt} must be < 1.\n'
,'Input 11: length(P)!=1.\n'
,'Input 12: length(delt)!=1.\n'
,'Input 13: length(alpha)!=1.\n'
,'Input 14: length(Xt) must be > 1.\n'
,'Input 15: Argument {Xs} must be of length 1.\n'
,'Input 16: Argument {Xt} must be of type vector!.\n'
,'Input 17: Argument {time} must be of type vector!.\n'
,'Input 18: Argument {alpha} must be of length 1.\n'
,'Input 19: Argument {Trunc} must be of length 2.\n'
,'Input 20: Cumulant truncation (Trunc[1]) must be either 4 or 8.\n'
,'Input 21: Density truncation (Trunc[2]) must be either 4 or 8.\n'
,'Input 22: Starting time {s} cannot be greater than {t}.\n'
,'Input 23: {Jdist} must be one of Normal, Exponential, Gamma or Laplace.\n'
,'Input 24: {Jtype} has to be of type Add or Mult.\n'
,'Input 25: {factorize} has to be TRUE or FALSE.\n'
,'Input 26: {factor.type} has to be Diffusion or Hawkes.\n'
,'Input 27: Diffusion terms are 0. Use {factorize=TRUE} and {factor.type="Hawkes"}.\n'
,'Input 28: Diffusion terms are 0. Density approximant cab perform poorly without \nfactorization. Try {factorize=TRUE}.\n'
)
warntrue=rep(FALSE,40)
if(length(Xt)<2){warntrue[14] =TRUE}#{stop(paste0(b1,warn[14],b2))}
if(length(Xs)!=1){warntrue[15] =TRUE}#{stop(paste0(b1,warn[15],b2))}
if(!is.vector(Xt)){warntrue[16] =TRUE}#{stop(paste0(b1,warn[16],b2))}
if(sum(Dindex)==0){warntrue[7] =TRUE}#{stop(paste0(b1,warn[7],b2))}
if(length(lower)>1){warntrue[1] =TRUE}#{stop(paste0(b1,warn[1],b2))}
if(length(upper)>1){warntrue[2] =TRUE}#{stop(paste0(b1,warn[2],b2))}
if(upper<=lower) {warntrue[3] =TRUE}# {stop(paste0(b1,warn[3],b2))}
#if((Dindex==3)|(Dindex==4)){if(lower[1]<=0){warntrue[4] =TRUE}}#{stop(paste0(b1,warn[4],b2))}}
if(Dindex==5){if(any(Xt<=0)|any(Xt>=1)){warntrue[5] =TRUE}}#{stop(paste0(b1,warn[5],b2))}}
if(P<10){warntrue[6] =TRUE}#{stop(paste0(b1,warn[6],b2))}
if(length(P)!=1){warntrue[11] =TRUE}#{stop(paste0(b1,warn[11],b2))}
if(length(alpha)!=1){warntrue[18] =TRUE}#{stop(paste0(b1,warn[18],b2))}
if(Trunc[2]>Trunc[1]){warntrue[8] =TRUE}#{stop(paste0(b1,warn[8],b2))}
if(length(Trunc)!=2){warntrue[19] =TRUE}#{stop(paste0(b1,warn[19],b2))}
if(sum(TR.order==c(4,8))!=1){warntrue[20] =TRUE}#{stop(paste0(b1,warn[20],b2))}
if(sum(DTR.order==c(4,8))!=1){warntrue[21] =TRUE}#{stop(paste0(b1,warn[21],b2))}
#if(delt>=1/5){stop(paste0(b1,warn[9],b2))}
if(delt>=1){warntrue[10] =TRUE}#{stop(paste0(b1,warn[10],b2))}
if(t<s){warntrue[22] =TRUE}#{stop(paste0(b1,warn[22],b2))}
if((Jtype!='Add')&&(Jtype!='Mult')) {warntrue[24] =TRUE}
if((factorize!=TRUE)&&(factorize!=FALSE)) {warntrue[25] =TRUE}
if((factor.type!='Diffusion')&&(factor.type!='Hawkes')) {warntrue[26] =TRUE}
if(sum(JDindex)==0){stop(paste0(b1,warn[23],b2))}
# Print output:
if(any(warntrue))
{
prnt = b1
for(i in which(warntrue))
{
prnt = paste0(prnt,warn[i])
}
prnt = paste0(prnt,b2)
stop(prnt)
}
diffusion.function=function()
{
time =seq(s,t,delt)[-1]
bod='0*x^0'
coefs=c('Q0','Q1','Q2')
obs=objects(pos=1)
for(i in 1:length(coefs))
{
if(sum(obs==coefs[i])){bod=paste0(bod,'+(',body(coefs[i])[2],')*x^',i-1)}
}
difff=function(x,t){}
body(difff) = parse(text=bod)
#print(body(difff))
M = matrix(0,length(Xt),length(time))
for(i in 1:length(time))
{
M[,i]=round(difff(Xt,time[i]),6)>0
}
if(sum(M)==length(time)*length(Xt)){return(1)}
if(sum(M)!=0) {return(2)}
if(sum(M)==0) {return(3)}
}
switch(diffusion.function(),
{},
{},
{
if((factorize)&&(factor.type!='Hawkes')){stop(paste0(b1,warn[27],b2))}
if(!factorize){warning(paste0(b1,warn[28],b2))}
})
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
function.list=objects(pos=1)
namess=c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta','Ja','Jb')
coef.terms2=paste0(namess,'(t)')
checklist=rep(0,length(namess))
for(i in 1:length(namess))
{
checklist[i]= sum(function.list== namess[i])
}
if(factorize)
{
if(sum(checklist[c(3,5,6)])==0){Dtype='Normal.A'}
}
if(Jdist=='Normal') {checklist[12:16]=0}
if(Jdist=='Exponential'){checklist[c(10,11,13,14,15,16)]=0}
if(Jdist=='Gamma') {checklist[c(10:12,15,16)]=0}
if(Jdist=='Laplace') {checklist[c(10:14)]=0}
if(Jdist=='Uniform') {checklist[c(10:14)]=0}
coef.index=which(checklist==1)
#==============================================================================
# Interface Module
#==============================================================================
namess4=namess
trim <- function (x) gsub("([[:space:]])", "", x)
for(i in 1:length(namess4))
{
if(sum(function.list==namess4[i]))
{
namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
}
}
namess4=matrix(namess4,length(namess4),1)
dinfo = c('Density approx. : ',
'Trunc. Order : ',
'Dens. Order : ')
dinfo[1] =paste0(dinfo[1],Dtype)
#if(Dindex!=1)
#{
# dinfo[2] =paste0(dinfo[2],P)
# dinfo[3] =paste0(dinfo[3],alpha)
#}
dinfo[2] =paste0(dinfo[2],Trunc[1])
dinfo[3] =paste0(dinfo[3],Trunc[2])
buffer0=c('================================================================')
buffer1=c('----------------------------------------------------------------')
buffer2=c('................................................................')
buffer3=c('... ... ... ... ... ... ... ... ... ... ... ')
buffer4=c('_____________________ Drift Coefficients _______________________')
buffer5=c('___________________ Diffusion Coefficients _____________________')
buffer6=c('__________________ Distribution Approximant ____________________')
buffer7=c('_______________________ Jump Mechanism _________________________')
buffer8=c('......................... Intensity ............................')
buffer9=c('........................... Jumps ..............................')
type.sol =" Jump Generalized Quadratic Diffusion (JGQD) "
Info=c(buffer0,type.sol,buffer0,buffer4,namess4[1:3],buffer5,namess4[4:6],buffer7,buffer8,namess4[7:9],buffer9,Jdist,namess4[9+which(checklist[-(1:9)]==1)],buffer6,dinfo)
Info=data.frame(matrix(Info,length(Info),1))
colnames(Info)=''
if(print.output)
{
print(Info,row.names = FALSE,right=F)
}
if(Integ)
{
if(TR.order ==8)
{
MAT=rbind(
c('1','1*x[1]','1*x[2]','','',''),
c('2*x[1]','2*x[2]','2*x[3]','1','1*x[1]','1*x[2]'),
c('3*x[2]','3*x[3]','3*x[4]','3*x[1]','3*x[2]','3*x[3]'),
c('4*x[3]','4*x[4]','4*x[5]','6*x[2]','6*x[3]','6*x[4]'),
c('5*x[4]','5*x[5]','5*x[6]','10*x[3]','10*x[4]','10*x[5]'),
c('6*x[5]','6*x[6]','6*x[7]','15*x[4]','15*x[5]','15*x[6]'),
c('7*x[6]','7*x[7]','7*x[8]','21*x[5]','21*x[6]','21*x[7]'),
c('8*x[7]','8*x[8]','8*(-(-9*x[1]*x[8]-36*x[2]*x[7]+72*x[1]^2*x[7]-84*x[3]*x[6]+504*x[1]*x[2]*x[6]-504*x[1]^3*x[6]-126*x[4]*x[5]+1008*x[1]*x[3]*x[5]+756*
x[2]^2*x[5]-4536*x[1]^2*x[2]*x[5]+3024*x[1]^4*x[5]+630*x[1]*x[4]^2+2520*x[2]*x[3]*x[4]-7560*x[1]^2*x[3]*x[4]-11340*x[1]*x[2]^2*x[4]+30240*x[1]^3*x[2]*x[4]
-15120*x[1]^5*x[4]+560*x[3]^3-15120*x[1]*x[2]*x[3]^2+20160*x[1]^3*x[3]^2-7560*x[2]^3*x[3]+90720*x[1]^2*x[2]^2*x[3]-151200*x[1]^4*x[2]*x[3]+60480*
x[1]^6*x[3]+22680*x[1]*x[2]^4-151200*x[1]^3*x[2]^3+272160*x[1]^5*x[2]^2-181440*x[1]^7*x[2]+40320*x[1]^9))','28*x[6]','28*x[7]','28*x[8]'))
MAT2=rbind(
c('1','1*x[8+1]','1*x[8+2]','','',''),
c('2*x[8+1]','2*x[8+2]','2*x[8+3]','1','1*x[8+1]','1*x[8+2]'),
c('3*x[8+2]','3*x[8+3]','3*x[8+4]','3*x[8+1]','3*x[8+2]','3*x[8+3]'),
c('4*x[8+3]','4*x[8+4]','4*x[8+5]','6*x[8+2]','6*x[8+3]','6*x[8+4]'),
c('5*x[8+4]','5*x[8+5]','5*x[8+6]','10*x[8+3]','10*x[8+4]','10*x[8+5]'),
c('6*x[8+5]','6*x[8+6]','6*x[8+7]','15*x[8+4]','15*x[8+5]','15*x[8+6]'),
c('7*x[8+6]','7*x[8+7]','7*x[8+8]','21*x[8+5]','21*x[8+6]','21*x[8+7]'),
c('8*x[8+7]','8*x[8+8]','8*(-(-9*x[8+1]*x[8+8]-36*x[8+2]*x[8+7]+72*x[8+1]^2*x[8+7]-84*x[8+3]*x[8+6]+504*x[8+1]*x[8+2]*x[8+6]-504*x[8+1]^3*x[8+6]-126*x[8+4]*x[8+5]+1008*x[8+1]*x[8+3]*x[8+5]+756*
x[8+2]^2*x[8+5]-4536*x[8+1]^2*x[8+2]*x[8+5]+3024*x[8+1]^4*x[8+5]+630*x[8+1]*x[8+4]^2+2520*x[8+2]*x[8+3]*x[8+4]-7560*x[8+1]^2*x[8+3]*x[8+4]-11340*x[8+1]*x[8+2]^2*x[8+4]+30240*x[8+1]^3*x[8+2]*x[8+4]
-15120*x[8+1]^5*x[8+4]+560*x[8+3]^3-15120*x[8+1]*x[8+2]*x[8+3]^2+20160*x[8+1]^3*x[8+3]^2-7560*x[8+2]^3*x[8+3]+90720*x[8+1]^2*x[8+2]^2*x[8+3]-151200*x[8+1]^4*x[8+2]*x[8+3]+60480*
x[8+1]^6*x[8+3]+22680*x[8+1]*x[8+2]^4-151200*x[8+1]^3*x[8+2]^3+272160*x[8+1]^5*x[8+2]^2-181440*x[8+1]^7*x[8+2]+40320*x[8+1]^9))','28*x[8+6]','28*x[8+7]','28*x[8+8]'))
dims=rep('',8)
dims2=rep('',8)
for(i in which(checklist[1:6]==1))
{
for(j in 1:8)
{
if(MAT[j,i]!='')
{
dims[j] = paste0(dims[j],'+(',body(namess[i])[2],')*',MAT[j,i])
dims2[j] = paste0(dims2[j],'+(',body(namess[i])[2],')*',MAT2[j,i])
}
}
}
if(Jtype=='Add')
{
JA0 =
c("+1*mm1"
,"+2*mm1*x[1]+1*mm2"
,"+3*mm1*x[2]+3*mm2*x[1]+1*mm3"
,"+4*mm1*x[3]+6*mm2*x[2]+4*mm3*x[1]+1*mm4"
,"+5*mm1*x[4]+10*mm2*x[3]+10*mm3*x[2]+5*mm4*x[1]+1*mm5"
,"+6*mm1*x[5]+15*mm2*x[4]+20*mm3*x[3]+15*mm4*x[2]+6*mm5*x[1]+1*mm6"
,"+7*mm1*x[6]+21*mm2*x[5]+35*mm3*x[4]+35*mm4*x[3]+21*mm5*x[2]+7*mm6*x[1]+1*mm7"
,"+8*mm1*x[7]+28*mm2*x[6]+56*mm3*x[5]+70*mm4*x[4]+56*mm5*x[3]+28*mm6*x[2]+8*mm7*x[1]+1*mm8")
JA1 =
c("+1*mm1*x[1]"
,"+2*mm1*x[2]+1*mm2*x[1]"
,"+3*mm1*x[3]+3*mm2*x[2]+1*mm3*x[1]"
,"+4*mm1*x[4]+6*mm2*x[3]+4*mm3*x[2]+1*mm4*x[1]"
,"+5*mm1*x[5]+10*mm2*x[4]+10*mm3*x[3]+5*mm4*x[2]+1*mm5*x[1]"
,"+6*mm1*x[6]+15*mm2*x[5]+20*mm3*x[4]+15*mm4*x[3]+6*mm5*x[2]+1*mm6*x[1]"
,"+7*mm1*x[7]+21*mm2*x[6]+35*mm3*x[5]+35*mm4*x[4]+21*mm5*x[3]+7*mm6*x[2]+1*mm7*x[1]"
,"+8*mm1*x[8]+28*mm2*x[7]+56*mm3*x[6]+70*mm4*x[5]+56*mm5*x[4]+28*mm6*x[3]+8*mm7*x[2]+1*mm8*x[1]")
JA2 =
c("+1*mm1*x[2]"
,"+2*mm1*x[3]+1*mm2*x[2]"
,"+3*mm1*x[4]+3*mm2*x[3]+1*mm3*x[2]"
,"+4*mm1*x[5]+6*mm2*x[4]+4*mm3*x[3]+1*mm4*x[2]"
,"+5*mm1*x[6]+10*mm2*x[5]+10*mm3*x[4]+5*mm4*x[3]+1*mm5*x[2]"
,"+6*mm1*x[7]+15*mm2*x[6]+20*mm3*x[5]+15*mm4*x[4]+6*mm5*x[3]+1*mm6*x[2]"
,"+7*mm1*x[8]+21*mm2*x[7]+35*mm3*x[6]+35*mm4*x[5]+21*mm5*x[4]+7*mm6*x[3]+1*mm7*x[2]"
,"+8*mm1*((-(-9*x[8+1]*x[8+8]-36*x[8+2]*x[8+7]+72*x[8+1]^2*x[8+7]-84*x[8+3]*x[8+6]+504*x[8+1]*x[8+2]*x[8+6]-504*x[8+1]^3*x[8+6]-126*x[8+4]*x[8+5]+1008*x[8+1]*x[8+3]*x[8+5]+756*
x[8+2]^2*x[8+5]-4536*x[8+1]^2*x[8+2]*x[8+5]+3024*x[8+1]^4*x[8+5]+630*x[8+1]*x[8+4]^2+2520*x[8+2]*x[8+3]*x[8+4]-7560*x[8+1]^2*x[8+3]*x[8+4]-11340*x[8+1]*x[8+2]^2*x[8+4]+30240*x[8+1]^3*x[8+2]*x[8+4]
-15120*x[8+1]^5*x[8+4]+560*x[8+3]^3-15120*x[8+1]*x[8+2]*x[8+3]^2+20160*x[8+1]^3*x[8+3]^2-7560*x[8+2]^3*x[8+3]+90720*x[8+1]^2*x[8+2]^2*x[8+3]-151200*x[8+1]^4*x[8+2]*x[8+3]+60480*
x[8+1]^6*x[8+3]+22680*x[8+1]*x[8+2]^4-151200*x[8+1]^3*x[8+2]^3+272160*x[8+1]^5*x[8+2]^2-181440*x[8+1]^7*x[8+2]+40320*x[8+1]^9))) +28*mm2*x[8]+56*mm3*x[7]+70*mm4*x[6]+56*mm5*x[5]+28*mm6*x[4]+8*mm7*x[3]+1*mm8*x[2]")
if((!missing(beta))&&(length(beta)==2))
{
a..=beta[1]
b..=beta[2]
JA0 =
c("(b..*x[1]+a..)*mm1",
"(b..^2*x[2]+2*a..*b..*x[1]+a..^2)*mm2+(2*b..*x[2]+2*a..*x[1])*mm1",
"(b..^3*x[3]+3*a..*b..^2*x[2]+3*a..^2*b..*x[1]+a..^3)*mm3+(3*b..^2*x[3]+6*a..*b..*x[2]+3*a..^2*x[1])*mm2+(3*b..*x[3]+3*a..*x[2])*mm1",
"(b..^4*x[4]+4*a..*b..^3*x[3]+6*a..^2*b..^2*x[2]+4*a..^3*b..*x[1]+a..^4)*mm4+(4*b..^3*x[4]+12*a..*b..^2*x[3]+12*a..^2*b..*x[2]+4*a..^3*x[1])*mm3+(6*b..^2*x[4]+12*a..*b..*x[3]+6*a..^2*x[2])*mm2+(4*b..*x[4]+4*a..*x[3])*mm1",
"(b..^5*x[5]*mm5+5*a..*b..^4*x[4]*mm5+10*a..^2*b..^3*x[3]*mm5+10*a..^3*b..^2*x[2]*mm5+5*a..^4*b..*x[1]*mm5+a..^5*mm5+5*b..^4*x[5]*mm4+20*a..*b..^3*x[4]*mm4+30*a..^2*b..^2*x[3]*mm4+20*a..^3*b..*x[2]*mm4+5*a..^4*x[1]*mm4+10*b..^3*x[5]*mm3+30*a..*b..^2*x[4]*mm3+30*a..^2*b..*x[3]*mm3+10*a..^3*x[2]*mm3+10*b..^2*x[5]*mm2+20*a..*b..*x[4]*mm2+10*a..^2*x[3]*mm2+5*b..*x[5]*mm1+5*a..*x[4]*mm1)",
"(b..^6*x[6]*mm6+6*a..*b..^5*x[5]*mm6+15*a..^2*b..^4*x[4]*mm6+20*a..^3*b..^3*x[3]*mm6+15*a..^4*b..^2*x[2]*mm6+6*a..^5*b..*x[1]*mm6+a..^6*mm6+6*b..^5*x[6]*mm5+30*a..*b..^4*x[5]*mm5+60*a..^2*b..^3*x[4]*mm5+60*a..^3*b..^2*x[3]*mm5
+30*a..^4*b..*x[2]*mm5+6*a..^5*x[1]*mm5+15*b..^4*x[6]*mm4+60*a..*b..^3*x[5]*mm4+90*a..^2*b..^2*x[4]*mm4+60*a..^3*b..*x[3]*mm4+15*a..^4*x[2]*mm4+20*b..^3*x[6]*mm3+60*a..*b..^2*x[5]*mm3+60*a..^2*b..*x[4]*mm3+20*a..^3*x[3]*mm3+15*
b..^2*x[6]*mm2+30*a..*b..*x[5]*mm2+15*a..^2*x[4]*mm2+6*b..*x[6]*mm1+6*a..*x[5]*mm1)",
"(b..^7*x[7]*mm7+7*a..*b..^6*x[6]*mm7+21*a..^2*b..^5*x[5]*mm7+35*a..^3*b..^4*x[4]*mm7+35*a..^4*b..^3*x[3]*mm7+21*a..^5*b..^2*x[2]*mm7+7*a..^6*b..*x[1]*mm7+a..^7*mm7+7*b..^6*x[7]*mm6+42*a..*b..^5*x[6]*mm6+105*a..^2*b..^4*x[5]*
mm6+140*a..^3*b..^3*x[4]*mm6+105*a..^4*b..^2*x[3]*mm6+42*a..^5*b..*x[2]*mm6+7*a..^6*x[1]*mm6+21*b..^5*x[7]*mm5+105*a..*b..^4*x[6]*mm5+210*a..^2*b..^3*x[5]*mm5+210*a..^3*b..^2*x[4]*mm5+105*a..^4*b..*x[3]*mm5+21*a..^5*x[2]*mm5+
35*b..^4*x[7]*mm4+140*a..*b..^3*x[6]*mm4+210*a..^2*b..^2*x[5]*mm4+140*a..^3*b..*x[4]*mm4+35*a..^4*x[3]*mm4+35*b..^3*x[7]*mm3+105*a..*b..^2*x[6]*mm3+105*a..^2*b..*x[5]*mm3+35*a..^3*x[4]*mm3+21*b..^2*x[7]*mm2+42*a..*b..*x[6]*
mm2+21*a..^2*x[5]*mm2+7*b..*x[7]*mm1+7*a..*x[6]*mm1)",
"(b..^8*x[8]*mm8+8*a..*b..^7*x[7]*mm8+28*a..^2*b..^6*x[6]*mm8+56*a..^3*b..^5*x[5]*mm8+70*a..^4*b..^4*x[4]*mm8+56*a..^5*b..^3*x[3]*mm8+28*a..^6*b..^2*x[2]*mm8+8*a..^7*b..*x[1]*mm8+a..^8*mm8+8*b..^7*x[8]*mm7+56*a..*b..^6*x[7]*mm7
+168*a..^2*b..^5*x[6]*mm7+280*a..^3*b..^4*x[5]*mm7+280*a..^4*b..^3*x[4]*mm7+168*a..^5*b..^2*x[3]*mm7+56*a..^6*b..*x[2]*mm7+8*a..^7*x[1]*mm7+28*b..^6*x[8]*mm6+168*a..*b..^5*x[7]*mm6+420*a..^2*b..^4*x[6]*mm6+560*a..^3*b..^3*x[5]*mm6
+420*a..^4*b..^2*x[4]*mm6+168*a..^5*b..*x[3]*mm6+28*a..^6*x[2]*mm6+56*b..^5*x[8]*mm5+280*a..*b..^4*x[7]*mm5+560*a..^2*b..^3*x[6]*mm5+560*a..^3*b..^2*x[5]*mm5+280*a..^4*b..*x[4]*mm5+56*a..^5*x[3]*mm5+70*b..^4*x[8]*mm4+280*
a..*b..^3*x[7]*mm4+420*a..^2*b..^2*x[6]*mm4+280*a..^3*b..*x[5]*mm4+70*a..^4*x[4]*mm4+56*b..^3*x[8]*mm3+168*a..*b..^2*x[7]*mm3+168*a..^2*b..*x[6]*mm3+56*a..^3*x[5]*mm3+28*b..^2*x[8]*mm2+56*a..*b..*x[7]*mm2+28*a..^2*x[6]*mm2+8
*b..*x[8]*mm1+8*a..*x[7]*mm1)")
print('Invoked special jump structure.')
}
if(checklist[7]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[7])[2],')*(',JA0[i],')')
}
}
if(checklist[8]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[8])[2],')*(',JA1[i],')')
}
}
if(checklist[9]==1)
{
for(i in 1:8)
{
dims[i] = paste0(dims[i],'+(',body(namess[9])[2],')*(',JA2[i],')')
}
}
}
#c('G0','G1','G2','Q0','Q1','Q2','Lam0','Lam1','Lam2','Jmu','Jsig','Jlam','Jalpha','Jbeta')
if(checklist[7]==0){Lam0=function(t){0}}
if(checklist[8]==0){Lam1=function(t){0}}
if(checklist[9]==0){Lam2=function(t){0}}
if(checklist[10]==0){Jmu=function(t){0}}
if(checklist[11]==0){Jsig=function(t){0}}
if(checklist[12]==0){Jlam=function(t){0}}
if(checklist[13]==0){Jalpha=function(t){0}}
if(checklist[14]==0){Jbeta=function(t){0}}
if(checklist[15]==0){Ja=function(t){0}}
if(checklist[16]==0){Jb=function(t){0}}
odekernel=paste0('c(',dims[1],'\n,',dims[2],'\n,',dims[3],'\n,',dims[4],'\n,',dims[5],'\n,',dims[6],'\n,',dims[7],'\n,',dims[8],'\n,',
dims2[1],'\n,',dims2[2],'\n,',dims2[3],'\n,',dims2[4],'\n,',dims2[5],'\n,',dims2[6],'\n,',dims2[7],'\n,',dims2[8],'\n,',
body(namess[7])[2],'\n,',paste('(',body(namess[8])[2],')*x[9]'),'\n,',paste('(',body(namess[9])[2],')*x[10]'),')')
}
if(TR.order ==4)
{
MAT=rbind(
c('1','1*x[1]','1*x[2]','','',''),
c('2*x[1]','2*x[2]','2*x[3]','1','1*x[1]','1*x[2]'),
c('3*x[2]','3*x[3]','3*x[4]','3*x[1]','3*x[2]','3*x[3]'),
c('4*x[3]','4*x[4]','4*x[5]','6*x[2]','6*x[3]','6*x[4]'))
MAT2=rbind(
c('1','1*x[4+1]','1*x[4+2]','','',''),
c('2*x[4+1]','2*x[4+2]','2*x[4+3]','1','1*x[4+1]','1*x[4+2]'),
c('3*x[4+2]','3*x[4+3]','3*x[4+4]','3*x[4+1]','3*x[4+2]','3*x[4+3]'),
c('4*x[4+3]','4*x[4+4]','4*x[4+5]','6*x[4+2]','6*x[4+3]','6*x[4+4]'))
dims=rep('',4)
dims2=rep('',4)
for(i in which(checklist[1:6]==1))
{
for(j in 1:4)
{
if(MAT[j,i]!='')
{
dims[j] = paste0(dims[j],'+(',body(namess[i])[2],')*',MAT[j,i])
dims2[j] = paste0(dims2[j],'+(',body(namess[i])[2],')*',MAT2[j,i])
}
}
}
if(Jtype=='Add')
{
JA0 =
c("+1*mm1"
,"+2*mm1*x[1]+1*mm2"
,"+3*mm1*x[2]+3*mm2*x[1]+1*mm3"
,"+4*mm1*x[3]+6*mm2*x[2]+4*mm3*x[1]+1*mm4")
JA1 =
c("+1*mm1*x[1]"
,"+2*mm1*x[2]+1*mm2*x[1]"
,"+3*mm1*x[3]+3*mm2*x[2]+1*mm3*x[1]"
,"+4*mm1*x[4]+6*mm2*x[3]+4*mm3*x[2]+1*mm4*x[1]")
JA2 =
c("+1*mm1*x[2]"
,"+2*mm1*x[3]+1*mm2*x[2]"
,"+3*mm1*x[4]+3*mm2*x[3]+1*mm3*x[2]"
,"+4*mm1*x[5]+6*mm2*x[4]+4*mm3*x[3]+1*mm4*x[2]")
}
if((!missing(beta))&&(length(beta)==2))
{
a..=beta[1]
b..=beta[2]
JA0 =
c("b..*x[1]*mm1+a..*mm1"
,"b..^2*x[2]*mm2+2*a..*b..*x[1]*mm2+a..^2*mm2+2*b..*x[2]*mm1+2*a..*x[1]*mm1"
,"b..^3*x[3]*mm3+3*a..*b..^2*x[2]*mm3+3*a..^2*b..*x[1]*mm3+a..^3*mm3+3*b..^2*x[3]*mm2+6*a..*b..*x[2]*mm2+3*a..^2*x[1]*mm2+3*b..*x[3]*mm1+3*a..*x[2]*mm1"
,"b..^4*x[4]*mm4+4*a..*b..^3*x[3]*mm4+6*a..^2*b..^2*x[2]*mm4+4*a..^3*b..*x[1]*mm4+a..^4*mm4+4*b..^3*x[4]*mm3+12*a..*b..^2*x[3]*mm3+12*a..^2*b..*x[2]*mm3+4*a..^3*x[1]*mm3+6*b..^2*x[4]*mm2+12*a..*b..*x[3]*mm2+6*a..^2*x[2]*mm2+4*b..*x[4]*mm1+4*a..*x[3]*mm1")
print('Invoked special jump structure.')
}
if(checklist[7]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[7])[2],')*(',JA0[i],')')
}
}
if(checklist[8]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[8])[2],')*(',JA1[i],')')
}
}
if(checklist[9]==1)
{
for(i in 1:4)
{
dims[i] = paste0(dims[i],'+(',body(namess[9])[2],')*(',JA2[i],')')
}
}
if(checklist[7]==0)
{
Lam0=function(t){0}
}
if(checklist[8]==0)
{
Lam1=function(t){0}
}
if(checklist[9]==0)
{
Lam2=function(t){0}
}
odekernel=paste0('c(',dims[1],'\n,',dims[2],'\n,',dims[3],'\n,',dims[4],'\n,',
dims2[1],'\n,',dims2[2],'\n,',dims2[3],'\n,',dims2[4],'\n,',
body(namess[7])[2],'\n,',paste('(',body(namess[8])[2],')*x[5]'),'\n,',paste('(',body(namess[9])[2],')*x[6]'),')')
}
if(Jdist=='Normal')
{
prem =
'mm1 = mu
mm2 = (mu^2)+ (sig^2)
mm3 = (mu^3)+ 3* (mu^1)*(sig^2)
mm4 = (mu^4)+ 6* (mu^2)*(sig^2)+3*(sig^4)
mm5 = (mu^5)+ 10*(mu^3)*(sig^2)+15*(mu^1)*(sig^4)
mm6 = (mu^6)+ 15*(mu^4)*(sig^2)+45*(mu^2)*(sig^4)+15*(sig^6)
mm7 = (mu^7)+ 21*(mu^5)*(sig^2)+105*(mu^3)*(sig^4)+105*mu*(sig^6)
mm8 = (mu^8)+ 28*(mu^6)*(sig^2)+210*(mu^4)*(sig^4)+420*(mu^2)*(sig^6)+105*(sig^8)
'
prem=paste0('mu = ',body('Jmu')[2],'\n','sig = ', body('Jsig')[2],'\n',prem)
}
if(Jdist=='Exponential')
{
prem =
'mm1=1*(mua^1)
mm2=2*(mua^2)
mm3=6*(mua^3)
mm4=24*(mua^4)
mm5=120*(mua^5)
mm6=720*(mua^6)
mm7=5040*(mua^7)
mm8=40320*(mua^8)
'
prem=paste0('mua = ',body('Jlam')[2],'\n',prem)
}
if(Jdist=='Gamma')
{
prem =
'mm1=(alph*bet)
mm2=(alph+1)*alph*bet^2
mm3=(alph+2)*(alph+1)*alph*bet^3
mm4=(alph+3)*(alph+2)*(alph+1)*alph*bet^4
mm5=(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^5
mm6=(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^6
mm7=(alph+6)*(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^7
mm8=(alph+7)*(alph+6)*(alph+5)*(alph+4)*(alph+3)*(alph+2)*(alph+1)*alph*bet^8
'
prem=paste0('alph = ',body('Jalpha')[2],'\n',prem)
prem=paste0('bet = ',body('Jbeta')[2],'\n',prem)
}
if(Jdist=='Laplace')
{
#print('yeah')
prem =
'mm1=0.5*(+2*aaaa^1*bbbb^0)
mm2=0.5*(+2*aaaa^2*bbbb^0+4*aaaa^0*bbbb^2)
mm3=0.5*(+2*aaaa^3*bbbb^0+12*aaaa^1*bbbb^2)
mm4=0.5*(+2*aaaa^4*bbbb^0+24*aaaa^2*bbbb^2+48*aaaa^0*bbbb^4)
mm5=0.5*(+2*aaaa^5*bbbb^0+40*aaaa^3*bbbb^2+240*aaaa^1*bbbb^4)
mm6=0.5*(+2*aaaa^6*bbbb^0+60*aaaa^4*bbbb^2+720*aaaa^2*bbbb^4+1440*aaaa^0*bbbb^6)
mm7=0.5*(+2*aaaa^7*bbbb^0+84*aaaa^5*bbbb^2+1680*aaaa^3*bbbb^4+10080*aaaa^1*bbbb^6)
mm8=0.5*(+2*aaaa^8*bbbb^0+112*aaaa^6*bbbb^2+3360*aaaa^4*bbbb^4+40320*aaaa^2*bbbb^6+80640*aaaa^0*bbbb^8)
'
prem=paste0('aaaa = ',body('Ja')[2],'\n',prem)
prem=paste0('bbbb = ',body('Jb')[2],'\n',prem)
}
if(Jdist=='Uniform')
{
prem =
'
mm1=(bbbb+aaaa)/2
mm2=(aaaa^2+aaaa*bbbb +bbbb^2)/3
mm3=1/4*(aaaa+bbbb)*(aaaa^2+bbbb^2)
mm4=1/5*(aaaa^4+aaaa^3*bbbb+aaaa^2*bbbb^2+aaaa*bbbb^3+bbbb^4)
mm5=1/6*(aaaa^5+aaaa^4*bbbb+aaaa^3*bbbb^2+aaaa^2*bbbb^3+aaaa^1*bbbb^4+bbbb^5)
mm6=1/7*(aaaa^6+aaaa^5*bbbb^1+aaaa^4*bbbb^2+aaaa^3*bbbb^3+aaaa^2*bbbb^4+aaaa^1*bbbb^5+bbbb^6)
mm7=1/8*(aaaa^7+aaaa^6*bbbb^1+aaaa^5*bbbb^2+aaaa^4*bbbb^3+aaaa^3*bbbb^4+aaaa^2*bbbb^5+aaaa^1*bbbb^6+bbbb^7)
mm8=1/9*(aaaa^8+aaaa^7*bbbb^1+aaaa^6*bbbb^2+aaaa^5*bbbb^3+aaaa^4*bbbb^4+aaaa^3*bbbb^5+aaaa^2*bbbb^6+aaaa^1*bbbb^7+bbbb^8)
'
prem=paste0('aaaa = ',body('Ja')[2],'\n',prem)
prem=paste0('bbbb = ',body('Jb')[2],'\n',prem)
}
odekernel=paste0('{',prem,odekernel,'}')
write(odekernel,'res.txt')
ff <- function(x,t){}
body(ff) = (parse(text=odekernel))
t.alpha=
c(0.000000000000000000000000000000000000000000000000000000000000,
0.100000000000000000000000000000000000000000000000000000000000,
0.539357840802981787532485197881302436857273449701009015505500,
0.809036761204472681298727796821953655285910174551513523258250,
0.309036761204472681298727796821953655285910174551513523258250,
0.981074190219795268254879548310562080489056746118724882027805,
0.833333333333333333333333333333333333333333333333333333333333,
0.354017365856802376329264185948796742115824053807373968324184,
0.882527661964732346425501486979669075182867844268052119663791,
0.642615758240322548157075497020439535959501736363212695909875,
0.357384241759677451842924502979560464040498263636787304090125,
0.117472338035267653574498513020330924817132155731947880336209,
0.833333333333333333333333333333333333333333333333333333333333,
0.309036761204472681298727796821953655285910174551513523258250,
0.539357840802981787532485197881302436857273449701009015505500,
0.100000000000000000000000000000000000000000000000000000000000,
1.00000000000000000000000000000000000000000000000000000000000)
N=round((t-s)/delt)
ttt=seq(s,t,by=delt)
solver=function()
{
MM=matrix(0,2*TR.order+3,N)
MA= c(rep(Xs^(1:TR.order),2),0,0,0)
tt=s
MM[,1]=MA
ii=0
while(ii < N)
{
ii=ii+1
x0=MA
fx0=ff(x0,tt)
x1=x0+delt*(0.1*fx0)
fx1=ff(x1,tt+t.alpha[2]*delt)
x2=x0+delt*(-0.915176561375291*fx0 +1.45453440217827*fx1)
fx2=ff(x2,tt+t.alpha[3]*delt)
x3=x0+delt*( 0.202259190301118*fx0 +0.606777570903354*fx2)
fx3=ff(x3,tt+t.alpha[4]*delt)
x4=x0+delt*( 0.184024714708644*fx0 +0.197966831227192*fx2-0.0729547847313633*fx3)
fx4=ff(x4,tt+t.alpha[5]*delt)
x5=x0+delt*( 0.0879007340206681*fx0 +0.410459702520261*fx3+0.482713753678866*fx4)
fx5=ff(x5,tt+t.alpha[6]*delt)
x6=x0+delt*(0.085970050490246*fx0 +0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5)
fx6=ff(x6,tt+t.alpha[7]*delt)
x7=x0+delt*(0.120930449125334*fx0 +0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6)
fx7=ff(x7,tt+t.alpha[8]*delt)
x8=x0+delt*(0.110854379580391*fx0 -0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7)
fx8=ff(x8,tt+t.alpha[9]*delt)
x9=x0+delt*(0.112054414752879*fx0 -0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8)
fx9=ff(x9,tt+t.alpha[10]*delt)
x10=x0+delt*(0.113976783964186*fx0 -0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9)
fx10=ff(x10,tt+t.alpha[11]*delt)
x11=x0+delt*(0.0798314528280196*fx0 -0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10)
fx11=ff(x11,tt+t.alpha[12]*delt)
x12=x0+delt*(0.985115610164857*fx0 +0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11)
fx12=ff(x12,tt+t.alpha[13]*delt)
x13=x0+delt*(0.895080295771633*fx0 +0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12)
fx13=ff(x13,tt+t.alpha[14]*delt)
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13)
fx14=ff(x14,tt+t.alpha[15]*delt)
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14)
fx15=ff(x15,tt+t.alpha[16]*delt)
x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15)
fx16=ff(x16,tt+t.alpha[17]*delt)
MA=MA+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt
tt=tt+delt
MM[,ii]=MA
setTxtProgressBar(pb,ii," "," ")
}
return(MM)
}
pb <- txtProgressBar(1,N,1,style = 1,width = 65)
MM=solver()
close(pb)
if(eval.density)
{
DD=matrix(0,length(Xt),N)
if(Dtype=='Saddle.Gamma')
{
print('yeaahh')
dens1=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
dens2=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p0 = exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
p1 = 1-p0
VVV[1,] = (MM[1,]-p0*MM[1+TR.order,])/p1
VVV[2,] = (MM[2,]-p0*MM[2+TR.order,])/p1
VVV[3,] = (MM[3,]-p0*MM[3+TR.order,])/p1
VVV[4,] = (MM[4,]-p0*MM[4+TR.order,])/p1
if(!factorize)
{
for(i in 1:N)
{
DD[,i]=dens1(Xt,MM[1:4,i])
}
}
if(factorize)
{
if(factor.type=='Diffusion')
{
for(i in 1:N)
{
if(p0[i]>0.01){DD[,i]=dens1(Xt,MMM[1:4,i])*p0[i]+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
if(factor.type=='Hawkes')
{
if(p0[i]>0.01){DD[,i]=p0[i]/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
DD=list(density=DD,MSH=0)
}
}
if(Dtype=='Gamma.Saddle')
{
print('yeahh')
dens2=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
dens1=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p0 = exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
p1 = 1-p0
VVV[1,] = (MM[1,]-p0*MM[1+TR.order,])/p1
VVV[2,] = (MM[2,]-p0*MM[2+TR.order,])/p1
VVV[3,] = (MM[3,]-p0*MM[3+TR.order,])/p1
VVV[4,] = (MM[4,]-p0*MM[4+TR.order,])/p1
if(!factorize)
{
for(i in 1:N)
{
DD[,i]=dens1(Xt,MM[1:4,i])
}
}
if(factorize)
{
if(factor.type=='Diffusion')
{
for(i in 1:N)
{
if(p0[i]>0.01){DD[,i]=dens1(Xt,MMM[1:4,i])*p0[i]+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
if(factor.type=='Hawkes')
{
if(p0[i]>0.01){DD[,i]=p0[i]/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+dens2(Xt,VVV[1:4,i])*p1[i]}else{DD[,i]=dens1(Xt,MM[1:4,i])}
}
}
DD=list(density=DD,MSH=0)
}
}
if(Dtype=='Gamma.A')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=dnorm(Xt,MMM[1,i],sqrt(MMM[2,i]-MMM[1,i]^2))*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
theta = (m[2]-m[1]^2)/m[1]
kappa = m[1]/theta
nu = dgamma(xx,shape=kappa,scale=theta)
return(nu)
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Normal.A')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=dnorm(Xt,MMM[1,i],sqrt(MMM[2,i]-MMM[1,i]^2))*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])/diff(Xt)[1]*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
return(dnorm(xx,k[1],sqrt(k[2])))
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Saddlepoint')
{
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=saddlep(Xt,MMM[1:4,i])*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^2)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^3)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,]^4)/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])*Xt*(abs(Xt-MM[1+TR.order,i])<0.99*diff(Xt)[1])+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
saddlep=function(xx,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
p=1/3*(3*(k[4]/6)*k[2] - ((k[3]/2)^2))/((k[4]/6)^2)
q=1/27*(27*((k[4]/6)^2)*(k[1]-xx) - 9*(k[4]/6)*(k[3]/2)*k[2] + 2*((k[3]/2)^3))/((k[4]/6)^3)
chk=(q^2)/4 + (p^3)/27
th=-(k[3]/2)/(3*(k[4]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)
K=k[1]*th+(k[2]*th^2)/2+(k[3]*th^3)/6 +(k[4]*th^4)/24
K1=k[1]+(k[2]*th)+(k[3]*th^2)/2+(k[4]*th^3)/6
K2=k[2]+(k[3]*th)+(k[4]*th^2)/2
K3=k[3]+(k[4]*th)
K4=k[4]
return(1/sqrt(2*pi*(K2))*exp(K-th*K1))
}
for(i in 1:N)
{
DD[,i]=saddlep(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order saddlepoint not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order saddlepoint not supported at present!')
}
}
if(Dtype=='Edgeworth')
{
edgeworth=function(x,m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
v3=k[3]/(sqrt(k[2])^3)
v4=k[4]/(sqrt(k[2])^4)
z=(x-k[1])/sqrt(k[2])
h2=function(z){z^2-1}
h3=function(z){z^3-3*z}
h4=function(z){z^4-6*z^2+3}
h5=function(z){z^5-10*z^3+15*z}
return(pmax(dnorm(x,k[1],sqrt(k[2]))*(1+v3/6*h3(z)+v4/24*h4(z)),0))
}
if(factorize)
{
if(factor.type=='Diffusion')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=edgeworth(Xt,MMM[1:4,i])*exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])+edgeworth(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(factor.type=='Hawkes')
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
p=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
VVV[1,] = MM[1,]#-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-p*MM[1+TR.order,]^2)/(1-p)
VVV[3,] = (MM[3,]-p*MM[1+TR.order,]^3)/(1-p)
VVV[4,] = (MM[4,]-p*MM[1+TR.order,]^4)/(1-p)
INDDDD = rep(0,N)
indexs=function(m)
{
k =m*0
k[1]= m[1]
k[2]= m[2]-1*k[1]*m[1]
k[3]= m[3]-1*k[1]*m[2]-2*k[2]*m[1]
k[4]= m[4]-1*k[1]*m[3]-3*k[2]*m[2]-3*k[3]*m[1]
skew = k[3]/(k[2]^(1.5))
kurt = k[4]/(k[2]^2)
return((skew^2+1)/kurt)
}
for(i in 1:N)
{
if(exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])>0.01)
{
DD[,i]=exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i])*Xt*(abs(Xt-MM[1,i])<0.005)+edgeworth(Xt,VVV[1:4,i])*(1-exp(-MM[2*TR.order+1,i]-MM[2*TR.order+2,i]-MM[2*TR.order+3,i]))#+dnorm(Xt,VVV[1,i],sqrt(VVV[2,i]-VVV[1,i]^2))*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))#+saddlep(Xt,VVV[1:4,i])*(1-exp(-MM[17,i]-MM[18,i]-MM[19,i]))
}else
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
INDDDD[i] = indexs(MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
}
if(!factorize)
{
if(DTR.order==4)
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
for(i in 1:N)
{
DD[,i]=edgeworth(Xt,MM[1:4,i])
}
DD=list(density=DD,MSH=0)
}
}
if(DTR.order==6)
{
stop('Incorrect input: 6th order Edgeworth not supported at present!')
}
if(DTR.order==8)
{
stop('Incorrect input: 8th order Edgeworth not supported at present!')
}
}
if(DTR.order==4)
{
pearson=function(x,k,type='Normal')
{
dx=diff(x)[1]
u = k*0
u[1,]=k[1,]
u[2,]=k[2,]
u[3,]=k[3,]
u[4,]=k[4,]
det=1*u[2,]*u[4,]+u[1,]*u[3,]*u[2,]+u[2,]*u[1,]*u[3,]-1*u[3,]*u[3,]-u[2,]*u[2,]*u[2,]-u[1,]*u[1,]*u[4,]
b11=(u[2,]*u[4,]-u[3,]*u[3,])/det
b12=(u[2,]*u[3,]-u[1,]*u[4,])/det
b13=(u[1,]*u[3,]-u[2,]*u[2,])/det
b21=(u[3,]*u[2,]-u[1,]*u[4,])/det
b22=(1*u[4,]-u[2,]*u[2,])/det
b23=(u[2,]*u[1,]-1*u[3,])/det
b31=(u[1,]*u[3,]-u[2,]*u[2,])/det
b32=(u[1,]*u[2,]-1*u[3,])/det
b33=(1*u[2,]-u[1,]*u[1,])/det
#Ntype
if(type=='Normal')
{
V1=0
V2=1
V3=2*u[1,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
thetas1=-betas1/1
thetas2=-betas2/2
thetas3=-betas3/3
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum(exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3)*rho*diff(y)[1])
DD[,i] = exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3)/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Gtype
if(type=='Gamma')
{
V1=1
V2=2*u[1,]
V3=3*u[2,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alphas =1-betas1
thetas1=-betas2/1
thetas2=-betas3/2
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum(tau^(alphas[i]-1)*exp(thetas1[i]*tau+thetas2[i]*tau^2)*rho*diff(y)[1])
DD[,i] = (x^(alphas[i]-1)*exp(thetas1[i]*x+thetas2[i]*x^2))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# IGtype
if(type=='InvGamma')
{
V1=2*u[1,]
V2=3*u[2,]
V3=4*u[3,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alphas =-betas2
nu =betas1
thetas =-betas3
DD=matrix(0,length(x),length(u[1,]))
MSH =matrix(0,P,length(u[1,]))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u[1,]))
{
y=seq(-1/(2*(u[1,i]-lower))*(sqrt(exp(2*alpha)+4*(u[1,i]-lower)^2)-exp(alpha)),-1/(2*(u[1,i]-upper))*(sqrt(exp(2*alpha)+4*(u[1,i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u[1,i]
MSH[2:(P-1),i] = tau
dtemp = sum((tau^alphas[i])*exp(nu[i]/tau+thetas[i]*tau)*rho*diff(y)[1])
DD[,i] = (x^alphas[i])*exp(nu[i]/x+thetas[i]*x)/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Btype
if(type=='Beta')
{
V1=1-2*u[1,]
V2=2*u[1,]-3*u[2,]
V3=3*u[2,]-4*u[3,]
betas1 =b11*V1+b12*V2+b13*V3
betas2 =b21*V1+b22*V2+b23*V3
betas3 =b31*V1+b32*V2+b33*V3
alpha =1-betas1
nu =1+(betas1+betas2+betas3)
theta1 =betas3
DD=matrix(0,length(x),length(u[1,]))
for(i in 1:length(x))
{
DD[i,]=(x[i]^(alpha-1))*((1-x[i])^(nu-1))*exp(theta1*x[i])
}
for(i in 1:dim(u)[2])
{
DD[,i]=DD[,i]/sum(DD[,i]*dx)
}
return(list(density=DD,MSH=0))
}
}
}
if(DTR.order==8)
{
pearson=function(x,k,type='Normal')
{
dx=diff(x)[1]
u1= k[1,]
u2=k[2,]
u3=k[3,]
u4=k[4,]
u5=k[5,]
u6=k[6,]
u7=k[7,]
u8=k[8,]
det=-u1*(u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*
(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))-u5*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5)))+u2*(u1*
(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*
(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4)))-u3*(u1*
(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))-u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u3*
(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))+u2*
(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))+u4*
(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u4*(u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*
(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*
(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))-u5*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))
b11= (u2*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))-u3*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))+u4*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))-u5*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5)))/det
b12= (-u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))+u2*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u3*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u4*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5)))/det
b13= (u1*(u3*(u6*u8-u7*u7)-u4*(u5*u8-u6*u7)+u5*(u5*u7-u6*u6))-u2*(u2*(u6*u8-u7*u7)-u4*(u4*u8-u5*u7)+u5*(u4*u7-u5*u6))+u3*(u2*(u5*u8-u6*u7)-u3*(u4*u8-u5*u7)+u5*(u4*u6-u5*u5))-u4*(u2*(u5*u7-u6*u6)-u3*(u4*u7-u5*u6)+u4*(u4*u6-u5*u5)))/det
b14= (-u1*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u6*u6)+u5*(u4*u7-u5*u6))+u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u5*u6)+u5*(u3*u7-u5*u5))-u3*(u2*(u4*u8-u6*u6)-u3*(u3*u8-u5*u6)+u5*(u3*u6-u4*u5))+u4*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u5*u5)+u4*(u3*u6-u4*u5)))/det
b15= (u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b21= (-u1*(u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)+u6*(u5*u7-u6*u6))+u3*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))-u4*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u5*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5)))/det
b22= (-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u3*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u4*(u6*u8-u7*u7)-u5*(u5*u8-u6*u7)-u4*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u6*(u5*u7-u6*u6))/det
b23= (u2*(u1*(u6*u8-u7*u7)-u4*(u3*u8-u4*u7)+u5*(u3*u7-u4*u6))-u3*(u1*(u5*u8-u6*u7)-u3*(u3*u8-u4*u7)+u5*(u3*u6-u4*u5))-u3*(u6*u8-u7*u7)+u4*(u5*u8-u6*u7)+u4*(u1*(u5*u7-u6*u6)-u3*(u3*u7-u4*u6)+u4*(u3*u6-u4*u5))-u5*(u5*u7-u6*u6))/det
b24= (-u2*(u1*(u5*u8-u6*u7)-u4*(u2*u8-u4*u6)+u5*(u2*u7-u4*u5))+u3*(u1*(u4*u8-u6*u6)-u3*(u2*u8-u4*u6)+u5*(u2*u6-u4*u4))+u3*(u5*u8-u6*u7)-u4*(u4*u8-u6*u6)-u4*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u4*u5)+u4*(u2*u6-u4*u4))+u5*(u4*u7-u5*u6))/det
b25= (u2*(u1*(u5*u7-u6*u6)-u4*(u2*u7-u3*u6)+u5*(u2*u6-u3*u5))-u3*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u3*u6)+u5*(u2*u5-u3*u4))-u3*(u5*u7-u6*u6)+u4*(u4*u7-u5*u6)+u4*(u1*(u4*u6-u5*u5)-u3*(u2*u6-u3*u5)+u4*(u2*u5-u3*u4))-u5*(u4*u6-u5*u5))/det
b31= (u1*(u3*(u6*u8-u7*u7)-u5*(u4*u8-u5*u7)+u6*(u4*u7-u5*u6))-u2*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))+u4*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u5*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4)))/det
b32= (u1*(u2*(u6*u8-u7*u7)-u5*(u3*u8-u4*u7)+u6*(u3*u7-u4*u6))-u3*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)-u3*(u6*u8-u7*u7)+u5*(u4*u8-u5*u7)+u4*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u6*(u4*u7-u5*u6))/det
b33= (-u1*(u1*(u6*u8-u7*u7)-u4*(u3*u8-u4*u7)+u5*(u3*u7-u4*u6))+u3*(u1*(u4*u8-u5*u7)-u2*(u3*u8-u4*u7)+u5*(u3*u5-u4*u4))+u2*(u6*u8-u7*u7)-u4*(u4*u8-u5*u7)-u4*(u1*(u4*u7-u5*u6)-u2*(u3*u7-u4*u6)+u4*(u3*u5-u4*u4))+u5*(u4*u7-u5*u6))/det
b34= (u1*(u1*(u5*u8-u6*u7)-u4*(u2*u8-u4*u6)+u5*(u2*u7-u4*u5))-u3*(u1*(u3*u8-u5*u6)-u2*(u2*u8-u4*u6)+u5*(u2*u5-u3*u4))-u2*(u5*u8-u6*u7)+u4*(u3*u8-u5*u6)+u4*(u1*(u3*u7-u5*u5)-u2*(u2*u7-u4*u5)+u4*(u2*u5-u3*u4))-u5*(u3*u7-u5*u5))/det
b35= (-u1*(u1*(u5*u7-u6*u6)-u4*(u2*u7-u3*u6)+u5*(u2*u6-u3*u5))+u3*(u1*(u3*u7-u4*u6)-u2*(u2*u7-u3*u6)+(u2*u4-u3*u3)*u5)+u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)-u4*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u3*u5)+u4*(u2*u4-u3*u3))+u5*(u3*u6-u4*u5))/det
b41= (-u1*(u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)+u6*(u4*u6-u5*u5))+u2*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))-u3*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)+u5*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b42= (-u1*(u2*(u5*u8-u6*u7)-u4*(u3*u8-u4*u7)+u6*(u3*u6-u4*u5))+u2*(u2*(u4*u8-u5*u7)-u3*(u3*u8-u4*u7)+(u3*u5-u4*u4)*u6)+u3*(u5*u8-u6*u7)-u4*(u4*u8-u5*u7)-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4))+u6*(u4*u6-u5*u5))/det
b43= (u1*(u1*(u5*u8-u6*u7)-u3*(u3*u8-u4*u7)+u5*(u3*u6-u4*u5))-u2*(u1*(u4*u8-u5*u7)-u2*(u3*u8-u4*u7)+u5*(u3*u5-u4*u4))-u2*(u5*u8-u6*u7)+u3*(u4*u8-u5*u7)+u4*(u1*(u4*u6-u5*u5)-u2*(u3*u6-u4*u5)+u3*(u3*u5-u4*u4))-u5*(u4*u6-u5*u5))/det
b44= (-u1*(u1*(u4*u8-u6*u6)-u3*(u2*u8-u4*u6)+u5*(u2*u6-u4*u4))+u2*(u1*(u3*u8-u5*u6)-u2*(u2*u8-u4*u6)+u5*(u2*u5-u3*u4))+u2*(u4*u8-u6*u6)-u3*(u3*u8-u5*u6)-u4*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u4*u4)+u3*(u2*u5-u3*u4))+u5*(u3*u6-u4*u5))/det
b45= (u1*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u3*u6)+u5*(u2*u5-u3*u4))-u2*(u1*(u3*u7-u4*u6)-u2*(u2*u7-u3*u6)+(u2*u4-u3*u3)*u5)-u2*(u4*u7-u5*u6)+u3*(u3*u7-u4*u6)+u4*(u1*(u3*u5-u4*u4)-u2*(u2*u5-u3*u4)+u3*(u2*u4-u3*u3))-u5*(u3*u5-u4*u4))/det
b51= (u1*(u3*(u5*u7-u6*u6)-u4*(u4*u7-u5*u6)+u5*(u4*u6-u5*u5))-u2*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))+u3*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u4*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4)))/det
b52= (u1*(u2*(u5*u7-u6*u6)-u4*(u3*u7-u4*u6)+u5*(u3*u6-u4*u5))-u2*(u2*(u4*u7-u5*u6)-u3*(u3*u7-u4*u6)+u5*(u3*u5-u4*u4))-u3*(u5*u7-u6*u6)+u4*(u4*u7-u5*u6)+u3*(u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)+u4*(u3*u5-u4*u4))-u5*(u4*u6-u5*u5))/det
b53= (-u1*(u1*(u5*u7-u6*u6)-u3*(u3*u7-u4*u6)+u4*(u3*u6-u4*u5))+u2*(u1*(u4*u7-u5*u6)-u2*(u3*u7-u4*u6)+u4*(u3*u5-u4*u4))+u2*(u5*u7-u6*u6)-u3*(u4*u7-u5*u6)-u3*(u1*(u4*u6-u5*u5)-u2*(u3*u6-u4*u5)+u3*(u3*u5-u4*u4))+u4*(u4*u6-u5*u5))/det
b54= (u1*(u1*(u4*u7-u5*u6)-u3*(u2*u7-u4*u5)+u4*(u2*u6-u4*u4))-u2*(u1*(u3*u7-u5*u5)-u2*(u2*u7-u4*u5)+u4*(u2*u5-u3*u4))-u2*(u4*u7-u5*u6)+u3*(u3*u7-u5*u5)+u3*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u4*u4)+u3*(u2*u5-u3*u4))-u4*(u3*u6-u4*u5))/det
b55= (-u1*(u1*(u4*u6-u5*u5)-u3*(u2*u6-u3*u5)+u4*(u2*u5-u3*u4))+u2*(u1*(u3*u6-u4*u5)-u2*(u2*u6-u3*u5)+u4*(u2*u4-u3*u3))+u2*(u4*u6-u5*u5)-u3*(u3*u6-u4*u5)-u3*(u1*(u3*u5-u4*u4)-u2*(u2*u5-u3*u4)+u3*(u2*u4-u3*u3))+u4*(u3*u5-u4*u4))/det
#Ntype
if(type=='Normal')
{
V1=0
V2=1
V3=2*u1
V4=3*u2
V5=4*u3
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
thetas1=-betas1/1
thetas2=-betas2/2
thetas3=-betas3/3
thetas4=-betas4/4
thetas5=-betas5/5
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum(exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3+thetas4[i]*tau^4+thetas5[i]*tau^5)*rho*diff(y)[1])
DD[,i] = (exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3+thetas4[i]*x^4+thetas5[i]*x^5))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Gtype
if(type=='Gamma')
{
V1=1
V2=2*u1
V3=3*u2
V4=4*u3
V5=5*u4
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alphas =1-betas1
thetas1=-betas2/1
thetas2=-betas3/2
thetas3=-betas4/3
thetas4=-betas5/4
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum(tau^(alphas[i]-1)*exp(thetas1[i]*tau+thetas2[i]*tau^2+thetas3[i]*tau^3+thetas4[i]*tau^4)*rho*diff(y)[1])
DD[,i] = (x^(alphas[i]-1)*exp(thetas1[i]*x+thetas2[i]*x^2+thetas3[i]*x^3+thetas4[i]*x^4))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# IGtype
if(type=='InvGamma')
{
V1=2*u1
V2=3*u2
V3=4*u3
V4=5*u4
V5=6*u5
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alphas =-betas2
nu =betas1
theta1 =-betas3
theta2 =-betas4/2
theta3 =-betas5/3
DD=matrix(0,length(x),length(u1))
MSH =matrix(0,P,length(u1))
MSH[1,] = lower
MSH[P,] = upper
for(i in 1:length(u1))
{
y=seq(-1/(2*(u1[i]-lower))*(sqrt(exp(2*alpha)+4*(u1[i]-lower)^2)-exp(alpha)),-1/(2*(u1[i]-upper))*(sqrt(exp(2*alpha)+4*(u1[i]-upper)^2)-exp(alpha)),length=P)
y=y[-c(1,length(y))]
rho = exp(alpha)*(1+y^2)/(1-(y)^2)^2
tau = exp(alpha)*(y)/(1-(y)^2)+u1[i]
MSH[2:(P-1),i] = tau
dtemp = sum((tau^alphas[i])*exp(nu[i]/tau+theta1[i]*tau+theta2[i]*tau^2+theta3[i]*tau^3)*rho*diff(y)[1])
DD[,i] = ((x^alphas[i])*exp(nu[i]/x+theta1[i]*x+theta2[i]*x^2+theta3[i]*x^3))/dtemp
}
return(list(density=DD,MSH=MSH))
}
# Btype
if(type=='Beta')
{
V1=1-2*u1
V2=2*u1-3*u2
V3=3*u2-4*u3
V4=4*u3-5*u4
V5=5*u4-6*u5
betas1 =b11*V1+b12*V2+b13*V3+b14*V4 +b15*V5
betas2 =b21*V1+b22*V2+b23*V3+b24*V4 +b25*V5
betas3 =b31*V1+b32*V2+b33*V3+b34*V4 +b35*V5
betas4 =b41*V1+b42*V2+b43*V3+b44*V4 +b45*V5
betas5 =b51*V1+b52*V2+b53*V3+b54*V4 +b55*V5
alpha =1-betas1
nu =1+(betas1+betas2+betas3+betas4+betas5)
theta1 =betas3+betas4+betas5
theta2 =(betas4+betas5)/2
theta3 =betas5/3
DD=matrix(0,length(x),length(u1))
for(i in 1:length(x))
{
DD[i,]=(x[i]^(alpha-1))*((1-x[i])^(nu-1))*exp(theta1*x[i]+theta2*x[i]^2+theta3*x[i]^3)
}
for(i in 1:length(u1))
{
DD[,i]=DD[,i]/sum(DD[,i]*dx)
}
return(list(density=DD,MSH=0))
}
}
}
if(Dtype=='Normal')
{
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Normal')
#print(d1)
d2=pearson(Xt,VVV,'Normal')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MMM,'Normal')
}
}
if(Dtype=='Gamma')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must be strictly positive for Dtype = Gamma. ')}
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Gamma')
#print(d1)
d2=pearson(Xt,VVV,'Gamma')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MMM,'Gamma')
}
}
if(Dtype=='InvGamma')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must be strictly positive for Dtype = InvGamma. ')}
DD=pearson(Xt,MM,'InvGamma')
}
if(Dtype=='Beta')
{
if(min(Xt)<=0){stop('Incorrect input: Xt must in [0,1] Dtype = Beta. Try seq(delta,1-delta,delta) for delta>0.')}
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
if(TR.order==8)
{
VVV[5,] = (MM[5,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[5+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[6,] = (MM[6,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[6+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[7,] = (MM[7,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[7+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[8,] = (MM[8,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[8+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
}
if(factorize)
{
prb=exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])
d1=pearson(Xt,MMM,'Beta')
#print(d1)
d2=pearson(Xt,VVV,'Beta')
DD=d1$density
#print(dim(DD))
for(i in 1:dim(d1$density)[2])
{
DD[,i]=d1$density[,i]*prb[i]+d2$density[,i]*(1-prb[i])
}
DD=list(density=DD)
}
if(!factorize)
{
DD=pearson(Xt,MM,'Beta')
}
}
}
}
if(!eval.density)
{
DD=list(density=NULL,MSH=NULL);
MMM = MM[(TR.order+1):(2*TR.order),]
VVV = MM*0
VVV[1,] = (MM[1,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[1+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[2,] = (MM[2,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[2+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[3,] = (MM[3,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[3+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
VVV[4,] = (MM[4,]-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,])*MM[4+TR.order,])/(1-exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]))
}
# Fix this!!!
if(TR.order==8)
{
u=MM[1:8,]
k=u*0
k[1,]= u[1,]
k[2,]= u[2,]-1*k[1,]*u[1,]
k[3,]= u[3,]-1*k[1,]*u[2,]-2*k[2,]*u[1,]
k[4,]= u[4,]-1*k[1,]*u[3,]-3*k[2,]*u[2,]-3*k[3,]*u[1,]
k[5,]= u[5,]-1*k[1,]*u[4,]-4*k[2,]*u[3,]-6*k[3,]*u[2,]-4*k[4,]*u[1,]
k[6,]= u[6,]-1*k[1,]*u[5,]-5*k[2,]*u[4,]-10*k[3,]*u[3,]-10*k[4,]*u[2,]-5*k[5,]*u[1,]
k[7,]= u[7,]-1*k[1,]*u[6,]-6*k[2,]*u[5,]-15*k[3,]*u[4,]-20*k[4,]*u[3,]-15*k[5,]*u[2,]-6*k[6,]*u[1,]
k[8,]= u[8,]-1*k[1,]*u[7,]-7*k[2,]*u[6,]-21*k[3,]*u[5,]-35*k[4,]*u[4,]-35*k[5,]*u[3,]-21*k[6,]*u[2,]-7*k[7,]*u[1,]
}
if(TR.order==4)
{
u=MM[1:4,]
k=u*0
k[1,]= u[1,]
k[2,]= u[2,]-1*k[1,]*u[1,]
k[3,]= u[3,]-1*k[1,]*u[2,]-2*k[2,]*u[1,]
k[4,]= u[4,]-1*k[1,]*u[3,]-3*k[2,]*u[2,]-3*k[3,]*u[1,]
}
ret= list(density=DD$density,Xt=Xt,time=seq(s,t,delt)[-1],cumulants=k,moments=MM,mesh=DD$MSH,zero_jump_prob =exp(-MM[2*TR.order+1,]-MM[2*TR.order+2,]-MM[2*TR.order+3,]),excess_moments =VVV)
class(ret) = 'JGQD.density'
return(ret)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.