# R/BiGQD.density.R In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions

```BiGQD.density=function(Xs,Ys,Xt,Yt,s,t,delt=1/100,Dtype='Saddlepoint',print.output=TRUE,eval.density=TRUE)
{
check_for_model=function()
{
txt=''
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
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()
a00=function(t){10+sin(2*pi*t)}
a10=function(t){-1}
c10=function(t){1}

b00=function(t){10+cos(2*pi*t)}
b01=function(t){-1}
f01=function(t){2}

model=BiGQD.density(10,10,seq(5,15,1/5),seq(5,15,1/5),0,5,delt=1/100)
GQD.plot(model)
--------------------------------------------------------------------------------
'
check=T
}
if((sum(func.list)>0)&&(sum(func.list[-c(1:12)])==0))
{
txt='
--------------------------------------------------------------------------------
At least two diffusion coefficients have to be defined! Try for example:
--------------------------------------------------------------------------------
GQD.remove()
c00=function(t){2+sin(2*pi*t)}
f00=function(t){2+sin(2*pi*t)}
model=BiGQD.density(10,10,seq(5,15,1/5),seq(5,15,1/5),0,3,delt=1/100)
GQD.plot(model)
--------------------------------------------------------------------------------
'
check=T
}
return(list(check=check,txt=txt))
}
check_for=check_for_model()
if(check_for[[1]]){stop(check_for[[2]])}

# Warning Module
Dindex = which(Dtypes==Dtype)
b1 = '\n==============================================================================\n'
b2 = '==============================================================================\n'
warn=c(
'Input 1: Dtype has to be one of Saddlepoint or Edgeworth.\n'
,'Input 2: Argument {delt} must be < 1.\n'
,'Input 3: length(delt)!=1.\n'
,'Input 4: length(Xt)=length(Yt) must be > 1.\n'
,'Input 5: Arguments {Xs,Ys} must be of length 1.\n'
,'Input 6: Arguments {Xt,Yt} must be of type vector!.\n'
,'Input 7: Starting time {s} cannot be greater than {t}.\n'
,'Input 8: length(Xt)!=length(Yt).\n'
,'Input 9: Arguments {s,t} must be of length 1.\n'
)

warntrue=rep(FALSE,40)

if(length(Xt)<2){warntrue[4] =TRUE}#{stop(paste0(b1,warn[4],b2))}
if(length(Xs)!=1){warntrue[5] =TRUE}#{stop(paste0(b1,warn[5],b2))}
if(length(Yt)<2){warntrue[4] =TRUE}#{stop(paste0(b1,warn[4],b2))}
if(length(Ys)!=1){warntrue[5] =TRUE}#{stop(paste0(b1,warn[5],b2))}
if(!is.vector(Xt)){warntrue[6] =TRUE}#{stop(paste0(b1,warn[6],b2))}
if(!is.vector(Yt)){warntrue[6] =TRUE}#{stop(paste0(b1,warn[6],b2))}
if(sum(Dindex)==0){warntrue[1] =TRUE}#{stop(paste0(b1,warn[1],b2))}
if(delt>=1){warntrue[2] =TRUE}#{stop(paste0(b1,warn[2],b2))}
if(length(delt)>1){warntrue[3] =TRUE}#{stop(paste0(b1,warn[3],b2))}
if(t<s){warntrue[7] =TRUE}#{stop(paste0(b1,warn[7],b2))}
if(length(Xt)!=length(Yt)){warntrue[8] =TRUE}#{stop(paste0(b1,warn[8],b2))}
if(length(t)!=1){warntrue[9] =TRUE}#{stop(paste0(b1,warn[9],b2))}
if(length(s)!=1){warntrue[9] =TRUE}#{stop(paste0(b1,warn[9],b2))}

# Print output:
if(any(warntrue))
{
prnt = b1
for(i in which(warntrue))
{
prnt = paste0(prnt,warn[i])
}
prnt = paste0(prnt,b2)
stop(prnt)
}

nnn=length(Xt)
X1=matrix(outer(Xt,rep(1,nnn)),1,nnn*nnn,byrow=T)
X2=matrix(t(outer(Yt,rep(1,nnn))),1,nnn*nnn,byrow=T)

pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
MAT=rbind(
c("()","(+1*k10)","(+1*k10*k10+1*k20)","(+1*k01)","(+1*k01*k01+1*k02)","(+1*k11+1*k10*k01)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+2*k20)","(+2*k10*k20+2*k20*k10+2*k30)","(+2*k11)","(+2*k01*k11+2*k11*k01+2*k12)","(+2*k21+2*k10*k11+2*k20*k01)","()","()","()","()","()","()","()","(+1*k10)","(+1*k20+1*k10*k10)","(+1*k01)","(+1*k02+1*k01*k01)","(+1*k11+1*k01*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+3*k30)","(+3*k10*k30+6*k20*k20+3*k30*k10+3*k40)","(+3*k21)","(+3*k01*k21+6*k11*k11+3*k21*k01+3*k22)","(+3*k31+3*k10*k21+6*k20*k11+3*k30*k01)","()","()","()","()","()","()","()","(+3*k20)","(+3*k30+3*k10*k20+3*k20*k10)","(+3*k11)","(+3*k12+3*k01*k11+3*k11*k01)","(+3*k21+3*k01*k20+3*k11*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+4*k40)","(+4*k10*k40+12*k20*k30+12*k30*k20+4*k40*k10)","(+4*k31)","(+4*k01*k31+12*k11*k21+12*k21*k11+4*k31*k01)","(+4*k10*k31+12*k20*k21+12*k30*k11+4*k40*k01)","()","()","()","()","()","()","()","(+6*k30)","(+6*k40+6*k10*k30+12*k20*k20+6*k30*k10)","(+6*k21)","(+6*k22+6*k01*k21+12*k11*k11+6*k21*k01)","(+6*k31+6*k01*k30+12*k11*k20+6*k21*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+1*k10)","(+1*k10*k10+1*k20)","(+1*k01)","(+1*k01*k01+1*k02)","(+1*k11+1*k01*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+2*k11)","(+2*k10*k11+2*k11*k10+2*k21)","(+2*k02)","(+2*k01*k02+2*k02*k01+2*k03)","(+2*k12+2*k01*k11+2*k02*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+1*k10)","(+1*k20+1*k10*k10)","(+1*k01)","(+1*k02+1*k01*k01)","(+1*k11+1*k01*k10)"
),c("()","()","()","()","()","()","()","(+3*k12)","(+3*k10*k12+6*k11*k11+3*k12*k10+3*k22)","(+3*k03)","(+3*k01*k03+6*k02*k02+3*k03*k01+3*k04)","(+3*k13+3*k01*k12+6*k02*k11+3*k03*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+3*k11)","(+3*k21+3*k10*k11+3*k11*k10)","(+3*k02)","(+3*k03+3*k01*k02+3*k02*k01)","(+3*k12+3*k01*k11+3*k02*k10)"
),c("()","()","()","()","()","()","()","(+4*k13)","(+4*k10*k13+12*k11*k12+12*k12*k11+4*k13*k10)","(+4*k04)","(+4*k01*k04+12*k02*k03+12*k03*k02+4*k04*k01)","(+4*k01*k13+12*k02*k12+12*k03*k11+4*k04*k10)","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+6*k12)","(+6*k22+6*k10*k12+12*k11*k11+6*k12*k10)","(+6*k03)","(+6*k04+6*k01*k03+12*k02*k02+6*k03*k01)","(+6*k13+6*k01*k12+12*k02*k11+6*k03*k10)"
),c("()","(+1*k11)","(+1*k10*k11+1*k11*k10+1*k21)","(+1*k02)","(+1*k01*k02+1*k02*k01+1*k03)","(+1*k12+1*k10*k02+1*k11*k01)","()","(+1*k20)","(+1*k10*k20+1*k20*k10+1*k30)","(+1*k11)","(+1*k01*k11+1*k11*k01+1*k12)","(+1*k21+1*k01*k20+1*k11*k10)","()","()","()","()","()","()","()","(+0.5*k10)","(+0.5*k20+0.5*k10*k10)","(+0.5*k01)","(+0.5*k02+0.5*k01*k01)","(+0.5*k11+0.5*k01*k10)","()","(+0.5*k10)","(+0.5*k20+0.5*k10*k10)","(+0.5*k01)","(+0.5*k02+0.5*k01*k01)","(+0.5*k11+0.5*k01*k10)","()","()","()","()","()","()"
),c("()","(+1*k12)","(+1*k10*k12+2*k11*k11+1*k12*k10+1*k22)","(+1*k03)","(+1*k01*k03+2*k02*k02+1*k03*k01+1*k04)","(+1*k13+1*k10*k03+2*k11*k02+1*k12*k01)","()","(+2*k21)","(+2*k10*k21+2*k11*k20+2*k20*k11+2*k21*k10+2*k31)","(+2*k12)","(+2*k01*k12+2*k02*k11+2*k11*k02+2*k12*k01+2*k13)","(+2*k22+2*k01*k21+2*k02*k20+2*k11*k11+2*k12*k10)","()","()","()","()","()","()","()","(+1*k11)","(+1*k21+1*k10*k11+1*k11*k10)","(+1*k02)","(+1*k03+1*k01*k02+1*k02*k01)","(+1*k12+1*k01*k11+1*k02*k10)","()","(+1*k11)","(+1*k21+1*k10*k11+1*k11*k10)","(+1*k02)","(+1*k03+1*k01*k02+1*k02*k01)","(+1*k12+1*k01*k11+1*k02*k10)","()","(+1*k20)","(+1*k30+1*k10*k20+1*k20*k10)","(+1*k11)","(+1*k12+1*k01*k11+1*k11*k01)","(+1*k21+1*k01*k20+1*k11*k10)"
),c("()","(+2*k21)","(+2*k10*k21+2*k11*k20+2*k20*k11+2*k21*k10+2*k31)","(+2*k12)","(+2*k01*k12+2*k02*k11+2*k11*k02+2*k12*k01+2*k13)","(+2*k22+2*k10*k12+2*k11*k11+2*k20*k02+2*k21*k01)","()","(+1*k30)","(+1*k10*k30+2*k20*k20+1*k30*k10+1*k40)","(+1*k21)","(+1*k01*k21+2*k11*k11+1*k21*k01+1*k22)","(+1*k31+1*k01*k30+2*k11*k20+1*k21*k10)","()","(+1*k11)","(+1*k21+1*k10*k11+1*k11*k10)","(+1*k02)","(+1*k03+1*k01*k02+1*k02*k01)","(+1*k12+1*k01*k11+1*k02*k10)","()","(+1*k20)","(+1*k30+1*k10*k20+1*k20*k10)","(+1*k11)","(+1*k12+1*k01*k11+1*k11*k01)","(+1*k21+1*k01*k20+1*k11*k10)","()","(+1*k20)","(+1*k30+1*k10*k20+1*k20*k10)","(+1*k11)","(+1*k12+1*k01*k11+1*k11*k01)","(+1*k21+1*k01*k20+1*k11*k10)","()","()","()","()","()","()"
),c("()","(+2*k22)","(+2*k10*k22+4*k11*k21+2*k12*k20+2*k20*k12+4*k21*k11+2*k22*k10)","(+2*k13)","(+2*k01*k13+4*k02*k12+2*k03*k11+2*k11*k03+4*k12*k02+2*k13*k01)","(+2*k10*k13+4*k11*k12+2*k12*k11+2*k20*k03+4*k21*k02+2*k22*k01)","()","(+2*k31)","(+2*k10*k31+2*k11*k30+4*k20*k21+4*k21*k20+2*k30*k11+2*k31*k10)","(+2*k22)","(+2*k01*k22+2*k02*k21+4*k11*k12+4*k12*k11+2*k21*k02+2*k22*k01)","(+2*k01*k31+2*k02*k30+4*k11*k21+4*k12*k20+2*k21*k11+2*k22*k10)","()","(+1*k12)","(+1*k22+1*k10*k12+2*k11*k11+1*k12*k10)","(+1*k03)","(+1*k04+1*k01*k03+2*k02*k02+1*k03*k01)","(+1*k13+1*k01*k12+2*k02*k11+1*k03*k10)","()","(+2*k21)","(+2*k31+2*k10*k21+2*k11*k20+2*k20*k11+2*k21*k10)","(+2*k12)","(+2*k13+2*k01*k12+2*k02*k11+2*k11*k02+2*k12*k01)","(+2*k22+2*k01*k21+2*k02*k20+2*k11*k11+2*k12*k10)","()","(+2*k21)","(+2*k31+2*k10*k21+2*k11*k20+2*k20*k11+2*k21*k10)","(+2*k12)","(+2*k13+2*k01*k12+2*k02*k11+2*k11*k02+2*k12*k01)","(+2*k22+2*k01*k21+2*k02*k20+2*k11*k11+2*k12*k10)","()","(+1*k30)","(+1*k40+1*k10*k30+2*k20*k20+1*k30*k10)","(+1*k21)","(+1*k22+1*k01*k21+2*k11*k11+1*k21*k01)","(+1*k31+1*k01*k30+2*k11*k20+1*k21*k10)"
),c("()","(+1*k13)","(+1*k10*k13+3*k11*k12+3*k12*k11+1*k13*k10)","(+1*k04)","(+1*k01*k04+3*k02*k03+3*k03*k02+1*k04*k01)","(+1*k10*k04+3*k11*k03+3*k12*k02+1*k13*k01)","()","(+3*k22)","(+3*k10*k22+6*k11*k21+3*k12*k20+3*k20*k12+6*k21*k11+3*k22*k10)","(+3*k13)","(+3*k01*k13+6*k02*k12+3*k03*k11+3*k11*k03+6*k12*k02+3*k13*k01)","(+3*k01*k22+6*k02*k21+3*k03*k20+3*k11*k12+6*k12*k11+3*k13*k10)","()","()","()","()","()","()","()","(+1.5*k12)","(+1.5*k22+1.5*k10*k12+3*k11*k11+1.5*k12*k10)","(+1.5*k03)","(+1.5*k04+1.5*k01*k03+3*k02*k02+1.5*k03*k01)","(+1.5*k13+1.5*k01*k12+3*k02*k11+1.5*k03*k10)","()","(+1.5*k12)","(+1.5*k22+1.5*k10*k12+3*k11*k11+1.5*k12*k10)","(+1.5*k03)","(+1.5*k04+1.5*k01*k03+3*k02*k02+1.5*k03*k01)","(+1.5*k13+1.5*k01*k12+3*k02*k11+1.5*k03*k10)","()","(+3*k21)","(+3*k31+3*k10*k21+3*k11*k20+3*k20*k11+3*k21*k10)","(+3*k12)","(+3*k13+3*k01*k12+3*k02*k11+3*k11*k02+3*k12*k01)","(+3*k22+3*k01*k21+3*k02*k20+3*k11*k11+3*k12*k10)"
),c("()","(+3*k31)","(+3*k10*k31+3*k11*k30+6*k20*k21+6*k21*k20+3*k30*k11+3*k31*k10)","(+3*k22)","(+3*k01*k22+3*k02*k21+6*k11*k12+6*k12*k11+3*k21*k02+3*k22*k01)","(+3*k10*k22+3*k11*k21+6*k20*k12+6*k21*k11+3*k30*k02+3*k31*k01)","()","(+1*k40)","(+1*k10*k40+3*k20*k30+3*k30*k20+1*k40*k10)","(+1*k31)","(+1*k01*k31+3*k11*k21+3*k21*k11+1*k31*k01)","(+1*k01*k40+3*k11*k30+3*k21*k20+1*k31*k10)","()","(+3*k21)","(+3*k31+3*k10*k21+3*k11*k20+3*k20*k11+3*k21*k10)","(+3*k12)","(+3*k13+3*k01*k12+3*k02*k11+3*k11*k02+3*k12*k01)","(+3*k22+3*k01*k21+3*k02*k20+3*k11*k11+3*k12*k10)","()","(+1.5*k30)","(+1.5*k40+1.5*k10*k30+3*k20*k20+1.5*k30*k10)","(+1.5*k21)","(+1.5*k22+1.5*k01*k21+3*k11*k11+1.5*k21*k01)","(+1.5*k31+1.5*k01*k30+3*k11*k20+1.5*k21*k10)","()","(+1.5*k30)","(+1.5*k40+1.5*k10*k30+3*k20*k20+1.5*k30*k10)","(+1.5*k21)","(+1.5*k22+1.5*k01*k21+3*k11*k11+1.5*k21*k01)","(+1.5*k31+1.5*k01*k30+3*k11*k20+1.5*k21*k10)","()","()","()","()","()","()"
))

MAT[1,1] ='1'           # ONES(N2)
MAT[5,7] ='1'
MAT[2,13]='1'
MAT[9,19]='1'
MAT[9,25]='1'
MAT[6,31]='1'

#==============================================================================
#                           Solution TYPE Module
#==============================================================================

namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')

func.list=rep(0,36)
obs=objects(pos=1)
for(i in 1:36)
{
if(sum(obs==namess[i])){func.list[i]=1}
}

set1.1=c(3,14,15,20,21)                          # X Nonlinear
set1.2=c(11,28,29,34,35)                         # Y Nonlinear
set4=c(3,5,6,9,11,12,14:18,20:24,26:30,32:36)    # Nonlinear terms
set5=c(4:6,8,9,12,16:18,22:24,26,27,30,32,33,36) # Dependant terms

state3.0=(     any(func.list[set1.1]==1)     &     any(func.list[set1.2]==1)     &(sum(func.list[set5])==0))  # Any nonlinear in BOTH     and NO interaction term
state3.1=((sum(any(func.list[set1.1]==1))==0)&(    any(func.list[set1.2]==1)    )&(sum(func.list[set5])==0))  # Any nonlinear Y ,Normal X and NO interaction term
state3.2=(     any(func.list[set1.1]==1)     &(sum(any(func.list[set1.2]==1))==0)&(sum(func.list[set5])==0))  # Any nonlinear X ,Normal Y and NO interaction term
state4  =(any(func.list[set4]==1)&any(func.list[set5]==1))                                                    # Any nonlinear term AND interaction terms
state1  =!(state3.0|state3.1|state3.2|state4)
sol.state=c(state1,state3.0,state3.1,state3.2,state4)
sol.state=c('2nd Ord. Truncation, Bivariate Normal','4th Ord. Truncation, Independent Saddlepoints ','2nd & 4th Ord. Truncation (Indep)','4th & 2nd Ord. Trunc. (Indep)',"4th Ord. Truncation, Bivariate-Saddlepoint")[which(sol.state==T)]

#==============================================================================
#                   Generate TYPE of Solution
#==============================================================================

colnames(MAT)= namess
func.list=rep(0,36)
obs=objects(pos=1)
for(i in 1:36)
{
if(sum(obs==namess[i])){func.list[i]=1}
}
dims=rep('(',14)
for(i in 1:14)
{
for(j in which(func.list==1))
{
if(MAT[i,j]!='()')
{
dims[i]=paste0(dims[i],'+(',body(namess[j])[2],')*',MAT[i,j])
}
}
dims[i]=paste0(dims[i],')')
}
if(any(dims=='()')){dims[which(dims=='()')]='(0)'}

bod='c('
bod=paste0(bod,dims[1])
for(i in 2:14)
{
bod=paste0(bod,',',dims[i])
}
bod=paste0(bod,')')
strng=c("k10=x[1];k20=x[2];k30=x[3];k40=x[4];k01=x[5];k02=x[6];k03=x[7];k04=x[8];k11=x[9];k12=x[10];k21=x[11];k22=x[12];k13=x[13];k31=x[14];")
res =paste0(strng,bod)
f=function(x,t){}
body(f)= as.call(c(as.name("{"), parse(text = res)))

trim <- function (x) gsub("([[:space:]])", "", x)
namess4=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
indnames =rep(0,36)
for(i in 1:36)
{
if(sum(obs==namess4[i]))
{
indnames[i]=TRUE
namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
}
}
namess4=matrix(namess4,length(namess4),1)
buffer0=c('================================================================')
buffer1=c('----------------------------------------------------------------')
buffer2=c('................................................................')
buffer3=c('...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ... ')
buffer4=c('_____________________ Drift Coefficients _______________________')
buffer5=c('___________________ Diffusion Coefficients _____________________')

#Info1=data.frame(rbind(buffer0,matrix(names4[1:6,],6,1),buffer3,matrix(names4[1:6+6,],6,1),buffer3,matrix(names4[1:6+12,],6,1),buffer3,matrix(names4[1:6+18,],6,1),buffer3,matrix(names4[1:6+24,],6,1),buffer3,matrix(names4[1:6+30,],6,1),buffer2,prior.list,buffer),check.names=F)
#colnames(Info1)=type.sol
#buffer00=data.frame(buffer0,check.names=F)
#colnames(buffer00)=''
#print(buffer00 ,row.names = FALSE,right=F)
#print(Info1,row.names = FALSE,right=F)
# print(indnames)
Info=c(buffer0,type.sol,buffer0,buffer4,
namess4[1:6][which(indnames[1:6]==T)],
buffer3,
namess4[7:12][which(indnames[7:12]==T)],
buffer5,
namess4[13:18][which(indnames[13:18]==T)],
buffer3,
namess4[19:24][which(indnames[19:24]==T)],
buffer3,
namess4[25:30][which(indnames[25:30]==T)],
buffer3,
namess4[31:36][which(indnames[31:36]==T)])
Info=data.frame(matrix(Info,length(Info),1))
colnames(Info)=''
if(print.output)
{
print(Info,row.names = FALSE,right=F)
}
N.mesh=(t-s)/delt+1
pb <- txtProgressBar(1,2*N.mesh,1,style = 1,width = 65)
#print(body(f))
solve.ode=function(xx0,yy0,TT,TT0,delt)
{
M=c(xx0,rep(0,3),yy0,rep(0,3+6))
N.mesh=(TT-TT0)/delt+1
MMMM=matrix(0,14,N.mesh)
MMMM[,1]=M
ttt=TT0
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)

for(i in 2:N.mesh)
{

x0=M
fx0=f(x0,ttt)
x1=x0+delt*(0.1*fx0)
fx1=f(x1,ttt+t.alpha[2]*delt)
x2=x0+delt*(-0.915176561375291*fx0  +1.45453440217827*fx1)
fx2=f(x2,ttt+t.alpha[3]*delt)
x3=x0+delt*( 0.202259190301118*fx0  +0.606777570903354*fx2)
fx3=f(x3,ttt+t.alpha[4]*delt)
x4=x0+delt*( 0.184024714708644*fx0  +0.197966831227192*fx2-0.0729547847313633*fx3)
fx4=f(x4,ttt+t.alpha[5]*delt)
x5=x0+delt*( 0.0879007340206681*fx0 +0.410459702520261*fx3+0.482713753678866*fx4)
fx5=f(x5,ttt+t.alpha[6]*delt)
x6=x0+delt*(0.085970050490246*fx0   +0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5)
fx6=f(x6,ttt+t.alpha[7]*delt)
x7=x0+delt*(0.120930449125334*fx0   +0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6)
fx7=f(x7,ttt+t.alpha[8]*delt)
x8=x0+delt*(0.110854379580391*fx0   -0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7)
fx8=f(x8,ttt+t.alpha[9]*delt)
x9=x0+delt*(0.112054414752879*fx0   -0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8)
fx9=f(x9,ttt+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=f(x10,ttt+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=f(x11,ttt+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=f(x12,ttt+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=f(x13,ttt+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=f(x14,ttt+t.alpha[15]*delt)
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14)
fx15=f(x15,ttt+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=f(x16,ttt+t.alpha[17]*delt)

M=M+(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
ttt=ttt+delt
MMMM[,i]=M

if(sum(is.na(M))>0){break;stop("ODEsolver - NAs produced.");}
setTxtProgressBar(pb,i," "," ")
}
return(MMMM)
}

MM=solve.ode(Xs,Ys,t,s,delt)

if(eval.density)
{
{
n1=dim(MM)[1]
n2=dim(MM)[2]

kkk=0
DDD=array(0,dim=c(nnn,nnn,n2))
a=rep(0,nnn*nnn)
b=rep(0,nnn*nnn)

for(lll in 2:n2)
{
setTxtProgressBar(pb,lll+N.mesh," "," ")
kkk=kkk+1
x=MM[,lll]

k10=x[1]
k20=x[2]
k30=x[3]
k40=x[4]
k01=x[5]
k02=x[6]
k03=x[7]
k04=x[8]
k11=x[9]
k12=x[10]
k21=x[11]
k22=x[12]
k13=x[13]
k31=x[14]
tme=proc.time()
ffab=function(xm1,xm2,a=rep(0,nnn*nnn),b=rep(0,nnn*nnn))
{
XMAT1=xm1
XMAT2=xm2

k=0
abser=rep(0.1,nnn*nnn)
while((k<5000)&(sum(abser>0.005)>0))
{
# K= k10*a+k01*b+1/2*k20*a^2+1/2*k02*b^2+1/6*k30*a^3+1/6*k03*b^3+1/24*k40*a^4+1/24*k04*b^4+k11*a*b+1/2*k12*a*b^2+1/2*k21*a^2*b+1/6*a*b^3*k13+1/6*a^3*b*k31+1/4*a^2*b^2*k22
gg=k10+k20*a+1/2*k30*a^2+1/6*k40*a^3 +k11*b +1/2*k12*b^2+k21*a*b+1/6*b^3*k13+1/2*a^2*b*k31+1/2*a*b^2*k22-XMAT1
hh=k01+k02*b+1/2*k03*b^2+1/6*k04*b^3 +k11*a +k12*a*b+1/2*k21*a^2+1/2*a*b^2*k13+1/6*a^3*k31+1/2*a^2*b*k22-XMAT2
gg1=k20+k30*a+1/2*k40*a^2+k21*b+a*b*k31+1/2*b^2*k22
gg2=k11 +k12*b+k21*a+1/2*b^2*k13+1/2*a^2*k31+a*b*k22
hh1=k11 +k12*b+k21*a+1/2*b^2*k13+1/2*a^2*k31+a*b*k22
hh2=k02+k03*b+1/2*k04*b^2 +k12*a+a*b*k13+1/2*a^2*k22
DK2=gg1*hh2-gg2*hh1
ar=(hh*gg2-gg*hh2)/(gg1*hh2-gg2*hh1)
br=(hh*gg1-gg*hh1)/(gg1*hh2-gg2*hh1)
anew  =a+ar
bnew  =b-br
abser =((anew-a)^2+(bnew-b)^2)
itrue =(abser>0.001)
ifalse=(abser<=0.001)
#print(itrue)
a=anew#*ifalse+anew*itrue
b=bnew#*ifalse+bnew*itrue
k=k+1
}

return(list(a=a,b=b,k=k))
}

fff=function(a,b,xmat1,xmat2)
{
K= k10*a+k01*b+1/2*k20*a^2+1/2*k02*b^2+1/6*k30*a^3+1/6*k03*b^3+1/24*k40*a^4+1/24*k04*b^4+k11*a*b+1/2*k12*a*b^2+1/2*k21*a^2*b+1/6*a*b^3*k13+1/6*a^3*b*k31+1/4*a^2*b^2*k22
gg=k10+k20*a+1/2*k30*a^2+1/6*k40*a^3 +k11*b +1/2*k12*b^2+k21*a*b+1/6*b^3*k13+1/2*a^2*b*k31+1/2*a*b^2*k22-xmat1
hh=k01+k02*b+1/2*k03*b^2+1/6*k04*b^3 +k11*a +k12*a*b+1/2*k21*a^2+1/2*a*b^2*k13+1/6*a^3*k31+1/2*a^2*b*k22-xmat2
gg1=k20+k30*a+1/2*k40*a^2+k21*b+a*b*k31+1/2*b^2*k22
gg2=k11 +k12*b+k21*a+1/2*b^2*k13+1/2*a^2*k31+a*b*k22
hh1=k11 +k12*b+k21*a+1/2*b^2*k13+1/2*a^2*k31+a*b*k22
hh2=k02+k03*b+1/2*k04*b^2 +k12*a+a*b*k13+1/2*a^2*k22
DK2=gg1*hh2-gg2*hh1
return(exp(K-a*xmat1-b*xmat2)/(2*pi)/sqrt(abs(DK2)))
}
dett = (k02*k20-k11^2)
a    = (k20*(k10-X1)-k11*(k01-X2))/dett
b    = (-k11*(k10-X1)+k02*(k01-X2))/dett
cc=ffab(X1,X2,a,b)

DDD[,,lll]=t(matrix(fff(cc\$a,cc\$b,X1,X2),nnn,nnn,byrow=T))
tme=proc.time()-tme
#print(tme)
}
}
if(Dtype=='Edgeworth')
{
ffff=function(xx,yy,m)
{
X=rep(xx,length(yy))
Y=rep(yy,each=length(xx))

H=function(x,i)
{
switch(i+1,
{1},
{x},
{x^2-1},
{x^3-3*x},
{x^4-6*x^2+3},
{x^5-10*x^3+15*x},
{x^6-15*x^4+45*x^2-15})
}

k10 =m[1]
k20 =m[2]
k30 =m[3]
k40 =m[4]
k01 =m[5]
k02 =m[6]
k03 =m[7]
k04 =m[8]
k11 =m[9]
k21 =m[10]
k12 =m[11]
k22 =m[12]
k13 =m[13]
k31 =m[14]

x = (X-k10)/sqrt(k20)
y = (Y-k01)/sqrt(k02)
rho20 =1
rho30 =k30/(sqrt(k20)^3)
rho40 =k40/(sqrt(k20)^4)

rho02 =1
rho03 =k03/(sqrt(k02)^3)
rho04 =k04/(sqrt(k02)^4)

rho11 =k11/(sqrt(k20)*sqrt(k02))
rho12 =k12/(sqrt(k20)*sqrt(k02)^2)
rho21 =k21/(sqrt(k20)^2*sqrt(k02))
rho22 =k22/(sqrt(k20)^2*sqrt(k02)^2)
rho13 =k13/(sqrt(k20)^1*sqrt(k02)^3)
rho31 =k31/(sqrt(k20)^3*sqrt(k02)^1)

rho = k11/sqrt(k02*k20)
dmvnorm <- function(x,y)
{
1/(2*pi*sqrt(k20*k02)*sqrt(1 - rho^2))*exp(-((x-k10)^2/k20 -2*rho*(x-k10)*(y-k01)/sqrt(k20 * k02)+(y-k01)^2/k02)/(2*(1-rho^2)))
}

A=H(x,3)*rho30+3*H(x,2)*H(y,1)*rho21+3*H(x,1)*H(y,2)*rho12+H(y,3)*rho03
B=H(x,4)*rho40+4*H(x,3)*H(y,1)*rho31+6*H(x,2)*H(y,2)*rho22+4*H(x,1)*H(y,3)*rho13+H(y,4)*rho04
CC=H(x,6)*rho30^2+6*H(x,5)*H(y,1)*rho21*rho30+6*H(x,4)*H(y,2)*rho12*rho30+2*H(x,3)*H(y,3)*rho03*rho30+9*H(x,4)*H(y,2)*rho21^2+18*H(x,3)*H(y,3)*rho12*rho21+6*H(x,2)*H(y,4)*rho03*rho21+9*H(x,2)*H(y,4)*rho12^2+6*H(x,1)*H(y,5)*rho03*rho12+H(y,6)*rho03^2
val=pmax(dmvnorm(X,Y)*(1+1/6*A+1/24*B+1/72*CC),0)

return(val)
}

DDD=array(0,dim=c(nnn,nnn,dim(MM)[2]))
DDD[1,1,1] = 1
for(i in 2:dim(MM)[2])
{
DDD[,,i]=ffff(Xt,Yt,MM[,i])
setTxtProgressBar(pb,i+N.mesh," "," ")
}

}
close(pb)
DD1=matrix(0,nnn,dim(MM)[2])
for(i in 1:length(Xt))
{
p=1/3*(3*(MM[4,]/6)*MM[2,] - ((MM[3,]/2)^2))/((MM[4,]/6)^2)
q=1/27*(27*((MM[4,]/6)^2)*(MM[1,]-Xt[i]) - 9*(MM[4,]/6)*(MM[3,]/2)*MM[2,] + 2*((MM[3,]/2)^3))/((MM[4,]/6)^3)
chk=abs((q^2)/4 + (p^3)/27)
th=-(MM[3,]/2)/(3*(MM[4,]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)

K=MM[1,]*th+(MM[2,]*th^2)/2+(MM[3,]*th^3)/6 +(MM[4,]*th^4)/24
K1=MM[1,]+(MM[2,]*th)+(MM[3,]*th^2)/2+(MM[4,]*th^3)/6
K2=MM[2,]+(MM[3,]*th)+(MM[4,]*th^2)/2
K3=MM[3,]+(MM[4,]*th)
K4=MM[4,]
DD1[i,]=1/sqrt(2*pi*abs(K2))*exp(K-th*K1)
}
DD2=matrix(0,nnn,dim(MM)[2])
for(i in 1:length(Yt))
{
p=1/3*(3*(MM[8,]/6)*MM[6,] - ((MM[7,]/2)^2))/((MM[8,]/6)^2)
q=1/27*(27*((MM[8,]/6)^2)*(MM[5,]-Yt[i]) - 9*(MM[8,]/6)*(MM[7,]/2)*MM[6,] + 2*((MM[7,]/2)^3))/((MM[8,]/6)^3)
chk=abs((q^2)/4 + (p^3)/27)
th=-(MM[7,]/2)/(3*(MM[8,]/6))+(-q/2+sqrt(chk))^(1/3)-(q/2+sqrt(chk))^(1/3)

K=MM[5,]*th+(MM[6,]*th^2)/2+(MM[7,]*th^3)/6 +(MM[8,]*th^4)/24
K1=MM[5,]+(MM[6,]*th)+(MM[7,]*th^2)/2+(MM[8,]*th^3)/6
K2=MM[6,]+(MM[7,]*th)+(MM[8,]*th^2)/2
K3=MM[7,]+(MM[8,]*th)
K4=MM[8,]
DD2[i,]=1/sqrt(2*pi*abs(K2))*exp(K-th*K1)
}
}

if(!eval.density)
{
DDD=NULL;DD1=NULL;DD2=NULL;setTxtProgressBar(pb,2*N.mesh," "," ");
}

rownames(MM) =c("k10","k20","k30","k40","k01","k02","k03","k04","k11","k12","k21","k22","k13","k31")
ret=list(density=DDD,Xmarginal=DD1,Ymarginal=DD2,Xt=Xt,Yt=Yt,time=seq(s,t,by=delt),cumulants=MM)
class(ret) = 'BiGQD.density'
return(ret)
}
```

## Try the DiffusionRgqd package in your browser

Any scripts or data that you put into this service are public.

DiffusionRgqd documentation built on May 2, 2019, 3:26 a.m.