#===============================================================================
# This example simulates a bivariate time homogeneous diffusion and shows how
# to conduct inference using BiGQD.mcmc(). We fit two competing models and then
# use the output to select a winner.
#-------------------------------------------------------------------------------
data(SDEsim2)
data(SDEsim2)
attach(SDEsim2)
# Have a look at the time series:
plot(Xt~time,type='l',col='blue',ylim=c(2,10),main='Simulated Data',xlab='Time (t)',ylab='State',
axes=FALSE)
lines(Yt~time,col='red')
expr1=expression(dX[t]==2(Y[t]-X[t])*dt+0.3*sqrt(X[t]*Y[t])*dW[t])
expr2=expression(dY[t]==(5-Y[t])*dt+0.5*sqrt(Y[t])*dB[t])
text(50,9,expr1)
text(50,8.5,expr2)
axis(1,seq(0,100,5))
axis(1,seq(0,100,5/10),tcl=-0.2,labels=NA)
axis(2,seq(0,20,2))
axis(2,seq(0,20,2/10),tcl=-0.2,labels=NA)
#------------------------------------------------------------------------------
# Define the coefficients of a proposed model
#------------------------------------------------------------------------------
GQD.remove()
a00 <- function(t){theta[1]*theta[2]}
a10 <- function(t){-theta[1]}
c00 <- function(t){theta[3]*theta[3]}
b00 <- function(t){theta[4]}
b01 <- function(t){-theta[5]}
f00 <- function(t){theta[6]*theta[6]}
theta.start <- c(3,3,3,3,3,3)
prop.sds <- c(0.15,0.16,0.04,0.99,0.19,0.04)
updates <- 50000
X <- cbind(Xt,Yt)
# Define prior distributions:
priors=function(theta){dunif(theta[1],0,100)*dunif(theta[4],0,100)}
# Run the MCMC procedure
m1=BiGQD.mcmc(X,time,10,theta.start,prop.sds,updates)
#------------------------------------------------------------------------------
# Remove old coefficients and define the coefficients of a new model
#------------------------------------------------------------------------------
GQD.remove()
a10 <- function(t){-theta[1]}
a01 <- function(t){theta[1]*theta[2]}
c11 <- function(t){theta[3]*theta[3]}
b00 <- function(t){theta[4]*theta[5]}
b01 <- function(t){-theta[4]}
f01 <- function(t){theta[6]*theta[6]}
theta.start <- c(3,3,3,3,3,3)
prop.sds <- c(0.16,0.02,0.01,0.18,0.12,0.01)
# Define prior distributions:
priors=function(theta){dunif(theta[1],0,100)*dunif(theta[4],0,100)}
# Run the MCMC procedure
m2=BiGQD.mcmc(X,time,10,theta.start,prop.sds,updates)
# Compare estimates:
GQD.estimates(m1)
GQD.estimates(m2)
#------------------------------------------------------------------------------
# Compare the two models
#------------------------------------------------------------------------------
GQD.dic(list(m1,m2))
#===============================================================================
Run the code above in your browser using DataLab