R/GQD.passage.R

In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions

GQD.passage<-function(Xs,B,theta,t,delt)
{
  #utils::globalVariables(c('fpt','ieq','solverFPT'),'DiffusionRgqd')
  fpt      =function(res,res2,N,delt){}
  ieq      =function(res1,res2,N,delt){}
  solverFPT=function(Xs,Bt,theta,N,delt,tt){}
  rm(list =c('fpt','ieq','solverFPT'))
  #solverFPT=function(Xs,Bt,theta,N,delt,tt){}
  #library(Rcpp)
  #library(RcppArmadillo)
  txt1=
    '
  
  #include <RcppArmadillo.h>
  
  // [[Rcpp::depends("RcppArmadillo")]]
  
  arma::vec f(arma::vec a,arma::vec theta,double t)
{
 
  arma::vec atemp(4);
  
  atemp(0)=(+(theta[1])*1+(theta[2])*a(0)+(theta[3])*(a(1)+a(0)*a(0)));
  atemp(1)=(+(theta[2])*2*a(1)+(theta[3])*(2*a(2)+4*a(0)*a(1))+(theta[4])*1+(theta[5])*a(0)+(theta[6])*(a(1)+a(0)*a(0)));
  atemp(2)=(+(theta[2])*3*a(2)+(theta[3])*(3*a(3)+6*a(0)*a(2)+6*a(1)*a(1))+(theta[5])*3*a(1)+(theta[6])*(3*a(2)+6*a(0)*a(1)));
  atemp(3)=(+(theta[2])*4*a(3)+(theta[3])*(8*a(0)*a(3)+24*a(1)*a(2))+(theta[5])*6*a(2)+(theta[6])*(6*a(3)+12*a(0)*a(2)+12*a(1)*a(1)));
  return atemp;
}
  
  // [[Rcpp::export]]
  double saddle(arma::vec xx,double Bt)
{
  
  double val=exp(-0.5*log(2*3.141592653589793*xx(1))-0.5*(Bt-xx(0))*(Bt-xx(0))/xx(1));
  return(val);
}
  // [[Rcpp::export]]
  double saddle2(arma::vec xx,double Bt)
{
  
  double val=-exp(-0.5*log(2*3.141592653589793*xx(1))-0.5*(Bt-xx(0))*(Bt-xx(0))/xx(1))*(Bt-xx(0))/xx(1);
  return(val);
}
  
  // [[Rcpp::export]]
  double pcurr(arma::vec xx,double Bt,arma::vec theta)
{
  double val =(0.5*(theta[1]+theta[2]*Bt+theta[3]*Bt*Bt)-0.75*(theta[5]+2*theta[6]*Bt))*saddle(xx,Bt)-0.5*(theta[4]+theta[5]*Bt+theta[6]*Bt*Bt)*saddle2(xx,Bt);
  return(val);
}
  
  // [[Rcpp::export]]
  arma::vec solverFPT(double Xs,double Bt,arma::vec theta,int N,double delt,double tt)
{
  arma::vec fx0(4);
  arma::vec fx1(4);
  arma::vec fx2(4);
  arma::vec fx3(4);
  arma::vec fx4(4);
  arma::vec fx5(4);
  arma::vec fx6(4);
  arma::vec fx7(4);
  arma::vec fx8(4);
  arma::vec fx9(4);
  arma::vec fx10(4);
  arma::vec fx11(4);
  arma::vec fx12(4);
  arma::vec fx13(4);
  arma::vec fx14(4);
  arma::vec fx15(4);
  arma::vec fx16(4);
  arma::vec x0(4);
  arma::vec x1(4);
  arma::vec x2(4);
  arma::vec x3(4);
  arma::vec x4(4);
  arma::vec x5(4);
  arma::vec x6(4);
  arma::vec x7(4);
  arma::vec x8(4);
  arma::vec x9(4);
  arma::vec x10(4);
  arma::vec x11(4);
  arma::vec x12(4);
  arma::vec x13(4);
  arma::vec x14(4);
  arma::vec x15(4);
  arma::vec x16(4);
  x0.fill(0);
  x0(0)=Xs;
  arma::vec fxx0(4);
  arma::vec fxx1(4);
  arma::vec fxx2(4);
  arma::vec fxx3(4);
  arma::vec fxx4(4);
  arma::vec fxx5(4);
  arma::vec fxx6(4);
  arma::vec fxx7(4);
  arma::vec fxx8(4);
  arma::vec fxx9(4);
  arma::vec fxx10(4);
  arma::vec fxx11(4);
  arma::vec fxx12(4);
  arma::vec fxx13(4);
  arma::vec fxx14(4);
  arma::vec fxx15(4);
  arma::vec fxx16(4);
  arma::vec xx0(4);
  arma::vec xx1(4);
  arma::vec xx2(4);
  arma::vec xx3(4);
  arma::vec xx4(4);
  arma::vec xx5(4);
  arma::vec xx6(4);
  arma::vec xx7(4);
  arma::vec xx8(4);
  arma::vec xx9(4);
  arma::vec xx10(4);
  arma::vec xx11(4);
  arma::vec xx12(4);
  arma::vec xx13(4);
  arma::vec xx14(4);
  arma::vec xx15(4);
  arma::vec xx16(4);
  xx0.fill(0);
  xx0(0)=Bt;
  arma::vec res1(N);
  arma::vec res2(N);
  
  double valA;
  double valB;
  double smtemp;
  double d=tt;
  for (int i = 1; i < N; i++)
{
  
fx0=f(x0,theta,d);
  x1=x0+delt*(0.1*fx0);
  fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
  fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
  fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
  x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
  fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
  fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
  x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
  fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
  fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
  x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
  fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
  x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
  fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*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,theta,d+0.357384241759677451842924502979560464040498263636787304090125*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,theta,d+0.117472338035267653574498513020330924817132155731947880336209*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,theta,d+0.833333333333333333333333333333333333333333333333333333333333*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,theta,d+0.309036761204472681298727796821953655285910174551513523258250*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,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
  fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*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,theta,d+delt);

  x0=x0+(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;

  valA=pcurr(x0,Bt,theta);
  
  fxx0=f(xx0,theta,d);
  xx1=xx0+delt*(0.1*fxx0);
  fxx1=f(xx1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  xx2=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1);
  fxx2=f(xx2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  xx3=xx0+delt*(0.202259190301118*fxx0+0.606777570903354*fxx2);
  fxx3=f(xx3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
  xx4=xx0+delt*(0.184024714708644*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3);
  fxx4=f(xx4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  xx5=xx0+delt*(0.0879007340206681*fxx0+0.410459702520261*fxx3+0.482713753678866*fxx4);
  fxx5=f(xx5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
  xx6=xx0+delt*(0.085970050490246*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-0.0731856375070851*fxx5);
  fxx6=f(xx6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  xx7=xx0+delt*(0.120930449125334*fxx0+0.260124675758296*fxx4+0.0325402621549091*fxx5-0.0595780211817361*fxx6);
  fxx7=f(xx7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
  xx8=xx0+delt*(0.110854379580391*fxx0-0.0605761488255006*fxx5+0.321763705601778*fxx6+0.510485725608063*fxx7);
  fxx8=f(xx8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
  xx9=xx0+delt*(0.112054414752879*fxx0-0.144942775902866*fxx5-0.333269719096257*fxx6+0.49926922955688*fxx7+0.509504608929686*fxx8);
  fxx9=f(xx9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
  xx10=xx0+delt*(0.113976783964186*fxx0-0.0768813364203357*fxx5+0.239527360324391*fxx6+0.397774662368095*fxx7+0.0107558956873607*fxx8-0.327769124164019*fxx9);
  fxx10=f(xx10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
  xx11=xx0+delt*(0.0798314528280196*fxx0-0.0520329686800603*fxx5-0.0576954146168549*fxx6+0.194781915712104*fxx7+0.145384923188325*fxx8-0.0782942710351671*fxx9-0.114503299361099*fxx10);
  fxx11=f(xx11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
  xx12=xx0+delt*(0.985115610164857*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-1.37896486574844*fxx5-0.861164195027636*fxx6+5.78428813637537*fxx7+3.28807761985104*fxx8-2.38633905093136*fxx9-3.25479342483644*fxx10-2.16343541686423*fxx11);
  fxx12=f(xx12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  xx13=xx0+delt*(0.895080295771633*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3-0.851236239662008*fxx5+0.398320112318533*fxx6+3.63937263181036*fxx7+1.5482287703983*fxx8-2.12221714704054*fxx9-1.58350398545326*fxx10-1.71561608285936*fxx11-0.0244036405750127*fxx12);
  fxx13=f(xx13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  xx14=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1+0*fxx2+0*fxx3-0.777333643644968*fxx4+0*fxx5-0.0910895662155176*fxx6+0.0910895662155176*fxx12+0.777333643644968*fxx13);
  fxx14=f(xx14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  xx15=xx0+delt*(0.1*fxx0-0.157178665799771*fxx2+0.157178665799771*fxx14);
  fxx15=f(xx15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  xx16=xx0+delt*(0.181781300700095*fxx0+0.675*fxx1+0.34275815984719*fxx2+0*fxx3+0.259111214548323*fxx4-0.358278966717952*fxx5-1.04594895940883*fxx6+0.930327845415627*fxx7+1.77950959431708*fxx8+0.1*fxx9-0.282547569539044*fxx10-0.159327350119973*fxx11-0.145515894647002*fxx12-0.259111214548323*fxx13-0.34275815984719*fxx14-0.675*fxx15);
  fxx16=f(xx16,theta,d+delt);

  xx0=xx0+(0.0333333333333333333333333333333333333333333333333333333333333*fxx0
        +0.0250000000000000000000000000000000000000000000000000000000000*fxx1
        +0.0333333333333333333333333333333333333333333333333333333333333*fxx2
        +0.000000000000000000000000000000000000000000000000000000000000*fxx3
        +0.0500000000000000000000000000000000000000000000000000000000000*fxx4
        +0.000000000000000000000000000000000000000000000000000000000000*fxx5
        +0.0400000000000000000000000000000000000000000000000000000000000*fxx6
        +0.000000000000000000000000000000000000000000000000000000000000*fxx7
        +0.189237478148923490158306404106012326238162346948625830327194*fxx8
        +0.277429188517743176508360262560654340428504319718040836339472*fxx9
        +0.277429188517743176508360262560654340428504319718040836339472*fxx10
        +0.189237478148923490158306404106012326238162346948625830327194*fxx11
        -0.0400000000000000000000000000000000000000000000000000000000000*fxx12
        -0.0500000000000000000000000000000000000000000000000000000000000*fxx13
        -0.0333333333333333333333333333333333333333333333333333333333333*fxx14
        -0.0250000000000000000000000000000000000000000000000000000000000*fxx15
        +0.0333333333333333333333333333333333333333333333333333333333333*fxx16)*delt;

   d=d+delt;
  valB=pcurr(xx0,Bt,theta);
  
  res1(i)=valB;
  if(i>1)
{
  smtemp=0;
  for (int j = 1; j < i; j++)
{
  smtemp=smtemp+res2(j-1)*res1(i-j);
}
  res2(i) = (2.0*valA-2.0*smtemp*delt);
}
}
  return(res2);
}
  
  
  '
  
  txt2=
    '
  
  #include <RcppArmadillo.h>
  
  // [[Rcpp::depends("RcppArmadillo")]]
  
  arma::vec f(arma::vec a,arma::vec theta,double t)
{
  
  arma::vec atemp(4);
  
  atemp(0)=(+(theta[1])*1+(theta[2])*a(0)+(theta[3])*(a(1)+a(0)*a(0)));
  atemp(1)=(+(theta[2])*2*a(1)+(theta[3])*(2*a(2)+4*a(0)*a(1))+(theta[4])*1+(theta[5])*a(0)+(theta[6])*(a(1)+a(0)*a(0)));
  atemp(2)=(+(theta[2])*3*a(2)+(theta[3])*(3*a(3)+6*a(0)*a(2)+6*a(1)*a(1))+(theta[5])*3*a(1)+(theta[6])*(3*a(2)+6*a(0)*a(1)));
  atemp(3)=(+(theta[2])*4*a(3)+(theta[3])*(8*a(0)*a(3)+24*a(1)*a(2))+(theta[5])*6*a(2)+(theta[6])*(6*a(3)+12*a(0)*a(2)+12*a(1)*a(1)));
  return atemp;
}
  
  // [[Rcpp::export]]
  double saddle(arma::vec xx,double Bt)
{
  double p=(1.0/3.0) *(3*(xx(3)/6.0)*xx(1) - pow(xx(2)/2.0,2))/pow(xx(3)/6.0,2);
  double q=(1.0/27.0)*(27*pow(xx(3)/6.0,2)*(xx(0)-Bt) - 9*(xx(3)/6.0)*(xx(2)/2.0)*xx(1) + 2*pow(xx(2)/2.0,3))/pow(xx(3)/6.0,3);
  double chk=pow(q,2)/4.0 + pow(p,3)/27.0;
  double th=-(xx(2)/2.0)/(3*(xx(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
  
  double K =xx(0)*th+(xx(1)*th*th)/2.0+(xx(2)*th*th*th)/6.0 +(xx(3)*th*th*th*th)/24.0;
  double K1=xx(0)   +(xx(1)*th)       +(xx(2)*th*th)/2.0    +(xx(3)*th*th*th)/6.0;
  double K2=xx(1)   +(xx(2)*th)       +(xx(3)*th*th)/2.0;
  double val=exp(-0.5*log(2*3.141592653589793*K2)+(K-th*K1));
  return(val);
}
  // [[Rcpp::export]]
  double saddle2(arma::vec xx,double Bt)
{
  double p=(1.0/3.0) *(3*(xx(3)/6.0)*xx(1) - pow(xx(2)/2.0,2))/pow(xx(3)/6.0,2);
  double q=(1.0/27.0)*(27*pow(xx(3)/6.0,2)*(xx(0)-Bt) - 9*(xx(3)/6.0)*(xx(2)/2.0)*xx(1) + 2*pow(xx(2)/2.0,3))/pow(xx(3)/6.0,3);
  double chk=pow(q,2)/4.0 + pow(p,3)/27.0;
  double th=-(xx(2)/2.0)/(3*(xx(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
  double thdash=1.0/(xx(1)+th*xx(2)+0.5*th*th*xx(3));
  
  double K =xx(0)*th+(xx(1)*th*th)/2.0+(xx(2)*th*th*th)/6.0 +(xx(3)*th*th*th*th)/24.0;
  double K1=xx(0)   +(xx(1)*th)       +(xx(2)*th*th)/2.0    +(xx(3)*th*th*th)/6.0;
  double K2=xx(1)   +(xx(2)*th)       +(xx(3)*th*th)/2.0;
  double gg      = 1.0/sqrt(2*3.141592653589793*(K2));
  double ggdash = -3.141592653589793*pow(2*3.141592653589793*(K2),-3.0/2.0)*(xx(2)*thdash+xx(3)*thdash*th);
  double hh      = exp(K-th*K1);
  double hhdash = exp(K-th*K1)*(-th*thdash*xx(1)-th*th*thdash*xx(2)-0.5*(th*th*th)*thdash*xx(3));
  
  double val=(ggdash*hh+gg*hhdash);
  return(val);
}
  
  // [[Rcpp::export]]
  double pcurr(arma::vec xx,double Bt,arma::vec theta)
{
  double val =(0.5*(theta[1]+theta[2]*Bt+theta[3]*Bt*Bt)-0.75*(theta[5]+2*theta[6]*Bt))*saddle(xx,Bt)-0.5*(theta[4]+theta[5]*Bt+theta[6]*Bt*Bt)*saddle2(xx,Bt);
  return(val);
}
  
  // [[Rcpp::export]]
  arma::vec solverFPT(double Xs,double Bt,arma::vec theta,int N,double delt,double tt)
{
arma::vec fx0(4);
  arma::vec fx1(4);
  arma::vec fx2(4);
  arma::vec fx3(4);
  arma::vec fx4(4);
  arma::vec fx5(4);
  arma::vec fx6(4);
  arma::vec fx7(4);
  arma::vec fx8(4);
  arma::vec fx9(4);
  arma::vec fx10(4);
  arma::vec fx11(4);
  arma::vec fx12(4);
  arma::vec fx13(4);
  arma::vec fx14(4);
  arma::vec fx15(4);
  arma::vec fx16(4);
  arma::vec x0(4);
  arma::vec x1(4);
  arma::vec x2(4);
  arma::vec x3(4);
  arma::vec x4(4);
  arma::vec x5(4);
  arma::vec x6(4);
  arma::vec x7(4);
  arma::vec x8(4);
  arma::vec x9(4);
  arma::vec x10(4);
  arma::vec x11(4);
  arma::vec x12(4);
  arma::vec x13(4);
  arma::vec x14(4);
  arma::vec x15(4);
  arma::vec x16(4);
  x0.fill(0);
  x0(0)=Xs;
  arma::vec fxx0(4);
  arma::vec fxx1(4);
  arma::vec fxx2(4);
  arma::vec fxx3(4);
  arma::vec fxx4(4);
  arma::vec fxx5(4);
  arma::vec fxx6(4);
  arma::vec fxx7(4);
  arma::vec fxx8(4);
  arma::vec fxx9(4);
  arma::vec fxx10(4);
  arma::vec fxx11(4);
  arma::vec fxx12(4);
  arma::vec fxx13(4);
  arma::vec fxx14(4);
  arma::vec fxx15(4);
  arma::vec fxx16(4);
  arma::vec xx0(4);
  arma::vec xx1(4);
  arma::vec xx2(4);
  arma::vec xx3(4);
  arma::vec xx4(4);
  arma::vec xx5(4);
  arma::vec xx6(4);
  arma::vec xx7(4);
  arma::vec xx8(4);
  arma::vec xx9(4);
  arma::vec xx10(4);
  arma::vec xx11(4);
  arma::vec xx12(4);
  arma::vec xx13(4);
  arma::vec xx14(4);
  arma::vec xx15(4);
  arma::vec xx16(4);
  xx0.fill(0);
  xx0(0)=Bt;
  arma::vec res1(N);
  arma::vec res2(N);
  
  double valA;
  double valB;
  double smtemp;
  double d=tt;
  for (int i = 1; i < N; i++)
{
  
  fx0=f(x0,theta,d);
  x1=x0+delt*(0.1*fx0);
  fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
  fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
  fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
  x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
  fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
  fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
  x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
  fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
  fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
  x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
  fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
  x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
  fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*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,theta,d+0.357384241759677451842924502979560464040498263636787304090125*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,theta,d+0.117472338035267653574498513020330924817132155731947880336209*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,theta,d+0.833333333333333333333333333333333333333333333333333333333333*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,theta,d+0.309036761204472681298727796821953655285910174551513523258250*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,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
  fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*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,theta,d+delt);

  x0=x0+(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;

  valA=pcurr(x0,Bt,theta);

  fxx0=f(xx0,theta,d);
  xx1=xx0+delt*(0.1*fxx0);
  fxx1=f(xx1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  xx2=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1);
  fxx2=f(xx2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  xx3=xx0+delt*(0.202259190301118*fxx0+0.606777570903354*fxx2);
  fxx3=f(xx3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
  xx4=xx0+delt*(0.184024714708644*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3);
  fxx4=f(xx4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  xx5=xx0+delt*(0.0879007340206681*fxx0+0.410459702520261*fxx3+0.482713753678866*fxx4);
  fxx5=f(xx5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
  xx6=xx0+delt*(0.085970050490246*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-0.0731856375070851*fxx5);
  fxx6=f(xx6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  xx7=xx0+delt*(0.120930449125334*fxx0+0.260124675758296*fxx4+0.0325402621549091*fxx5-0.0595780211817361*fxx6);
  fxx7=f(xx7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
  xx8=xx0+delt*(0.110854379580391*fxx0-0.0605761488255006*fxx5+0.321763705601778*fxx6+0.510485725608063*fxx7);
  fxx8=f(xx8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
  xx9=xx0+delt*(0.112054414752879*fxx0-0.144942775902866*fxx5-0.333269719096257*fxx6+0.49926922955688*fxx7+0.509504608929686*fxx8);
  fxx9=f(xx9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
  xx10=xx0+delt*(0.113976783964186*fxx0-0.0768813364203357*fxx5+0.239527360324391*fxx6+0.397774662368095*fxx7+0.0107558956873607*fxx8-0.327769124164019*fxx9);
  fxx10=f(xx10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
  xx11=xx0+delt*(0.0798314528280196*fxx0-0.0520329686800603*fxx5-0.0576954146168549*fxx6+0.194781915712104*fxx7+0.145384923188325*fxx8-0.0782942710351671*fxx9-0.114503299361099*fxx10);
  fxx11=f(xx11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
  xx12=xx0+delt*(0.985115610164857*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-1.37896486574844*fxx5-0.861164195027636*fxx6+5.78428813637537*fxx7+3.28807761985104*fxx8-2.38633905093136*fxx9-3.25479342483644*fxx10-2.16343541686423*fxx11);
  fxx12=f(xx12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
  xx13=xx0+delt*(0.895080295771633*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3-0.851236239662008*fxx5+0.398320112318533*fxx6+3.63937263181036*fxx7+1.5482287703983*fxx8-2.12221714704054*fxx9-1.58350398545326*fxx10-1.71561608285936*fxx11-0.0244036405750127*fxx12);
  fxx13=f(xx13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
  xx14=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1+0*fxx2+0*fxx3-0.777333643644968*fxx4+0*fxx5-0.0910895662155176*fxx6+0.0910895662155176*fxx12+0.777333643644968*fxx13);
  fxx14=f(xx14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
  xx15=xx0+delt*(0.1*fxx0-0.157178665799771*fxx2+0.157178665799771*fxx14);
  fxx15=f(xx15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
  xx16=xx0+delt*(0.181781300700095*fxx0+0.675*fxx1+0.34275815984719*fxx2+0*fxx3+0.259111214548323*fxx4-0.358278966717952*fxx5-1.04594895940883*fxx6+0.930327845415627*fxx7+1.77950959431708*fxx8+0.1*fxx9-0.282547569539044*fxx10-0.159327350119973*fxx11-0.145515894647002*fxx12-0.259111214548323*fxx13-0.34275815984719*fxx14-0.675*fxx15);
  fxx16=f(xx16,theta,d+delt);

  xx0=xx0+(0.0333333333333333333333333333333333333333333333333333333333333*fxx0
        +0.0250000000000000000000000000000000000000000000000000000000000*fxx1
        +0.0333333333333333333333333333333333333333333333333333333333333*fxx2
        +0.000000000000000000000000000000000000000000000000000000000000*fxx3
        +0.0500000000000000000000000000000000000000000000000000000000000*fxx4
        +0.000000000000000000000000000000000000000000000000000000000000*fxx5
        +0.0400000000000000000000000000000000000000000000000000000000000*fxx6
        +0.000000000000000000000000000000000000000000000000000000000000*fxx7
        +0.189237478148923490158306404106012326238162346948625830327194*fxx8
        +0.277429188517743176508360262560654340428504319718040836339472*fxx9
        +0.277429188517743176508360262560654340428504319718040836339472*fxx10
        +0.189237478148923490158306404106012326238162346948625830327194*fxx11
        -0.0400000000000000000000000000000000000000000000000000000000000*fxx12
        -0.0500000000000000000000000000000000000000000000000000000000000*fxx13
        -0.0333333333333333333333333333333333333333333333333333333333333*fxx14
        -0.0250000000000000000000000000000000000000000000000000000000000*fxx15
        +0.0333333333333333333333333333333333333333333333333333333333333*fxx16)*delt;

   d=d+delt;
   valB=pcurr(xx0,Bt,theta);

  
  res1(i)=valB;
  if(i>1)
{
  smtemp=0;
  for (int j = 1; j < i; j++)
{
  smtemp=smtemp+res2(j-1)*res1(i-j);
}
  res2(i) = (2.0*valA-2.0*smtemp*delt);
}
}
  return(res2);
}
  
  
  '

     if(sum(theta[c(3,5,6)] ==0)==3)
     {
       stre="Compiling C++ code. Please wait."
       cat(stre, " \r")
       flush.console()
       sourceCpp(code=txt1)
       cat('                                     ','\r')
     }
     if(sum(theta[c(3,5,6)] ==0)!=3)
     {
       stre="Compiling C++ code. Please wait."
       cat(stre, " \r")
       flush.console()
       sourceCpp(code=txt2)
       cat('                                     ','\r')
       
     }
     N=round(t/delt)
     res=solverFPT(Xs,B,c(0,theta),N,delt,0)
     return(list(density=res, time = seq(0,t-delt,delt),prob.coverage = round(sum(res*delt),3)))
 }