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
  Dtypes =c('Saddlepoint','Edgeworth')
  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)))
  
  type.sol ="                  GENERALIZED QUADRATIC DIFFUSON"
  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)
 {
  if(Dtype=='Saddlepoint')
  {
    n1=dim(MM)[1]
    n2=dim(MM)[2]
    
    kkk=0
    DDD=array(0,dim=c(nnn,nnn,n2))
    nr.updates=rep(0,n2-1)
    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)

      nr.updates[lll-1]=cc$k
      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)
    }