# SDEexamples.ssc examples of simulating from various SDEs # # author: Eric Zivot # created: March 2, 2004 # updated: March 25, 2005 # # plot 3-month T-bill data plot(tbill.dat[,"tb3mo"],reference.grid=F) # Repeat this for seed = 0, 1, 2, 3; use tsplot for the first plot, # then lines for the subsequent plots. Finish with abline(). sim0 = euler.pcode.gensim(rho=c(.5, 5), n.sim = 1000, n.var = 1, aux=euler.pcode.aux(ndt = 25, t.per.sim = 1, N=1, M=1, X0=0, seed = 0, ubound = 1000, drift.expr=expression(rho[1]), diffuse.expr=expression(rho[2]))) tsplot(sim) sim1 = euler.pcode.gensim(rho=c(.5, 5), n.sim = 1000, n.var = 1, aux=euler.pcode.aux(ndt = 25, t.per.sim = 1, N=1, M=1, X0=0, seed = 1, ubound = 1000, drift.expr=expression(rho[1]), diffuse.expr=expression(rho[2]))) lines(sim1) sim2 = euler.pcode.gensim(rho=c(.5, 5), n.sim = 1000, n.var = 1, aux=euler.pcode.aux(ndt = 25, t.per.sim = 1, N=1, M=1, X0=0, seed = 2, ubound = 1000, drift.expr=expression(rho[1]), diffuse.expr=expression(rho[2]))) lines(sim2) sim3 = euler.pcode.gensim(rho=c(.5, 5), n.sim = 1000, n.var = 1, aux=euler.pcode.aux(ndt = 25, t.per.sim = 1, N=1, M=1, X0=0, seed = 3, ubound = 1000, drift.expr=expression(rho[1]), diffuse.expr=expression(rho[2]))) lines(sim3) abline(0, .5) # simulate from OU process # For OU.gensim, rho is assumed packed as c(kappa, theta, sigma) kappa = .4; theta = .08; sigma = .1 sim.ou = OU.gensim(rho = c(kappa, theta, sigma), n.sim = 1000, n.burn = 500, aux = OU.aux(X0 = .1, ndt = 25, seed = 1, t.per.sim = 1/52)) # Plot the simulation tsplot(sim.ou) # Plot the long-run mean abline(h=theta) # simulate from CIR process # rho is packed as c(kappa, theta, sigma) kappa = .4; theta = .08; sigma = .1 rho = c(kappa, theta, sigma) # For this example we will pass in random normals for Brownian motion, # although CIR.gensim can generate these internally on the fly. ndt = 25; n.sim = 1000; n.burn = 500 set.seed(1) z = rnorm(ndt*(n.burn + n.sim)) sim.cir = CIR.gensim(rho = rho, n.sim = n.sim, n.burn = n.burn, aux = CIR.aux(X0 = theta, ndt = ndt, z = z, t.per.sim = 1/52)) # Plot the simulation tsplot(sim.cir) # Plot the long-run mean abline(h=theta) # CIR model can generate explosive simulations! set.seed(10) sim.cir3 = euler.cir(kappa = .4, theta = .08, sigma = 1, n.sim = 1000, n.burn = 0, X0 = .1, ndt = 100, t.per.sim = 1/52) # Plot shows explosive segment. tsplot(sim.cir3) # This using Doss transform is terrible! set.seed(0); z = rnorm(20*1000) sim.cir.a = atsm.fac.sim(Kappa = 1, Theta = .4, Sigma = .1, alpha = 0, beta = 1, form = "can", N = 1, m = 1, nsim = 1000, n0 = 0, ndt = 20, t.per.sim = .01, X0 = .1, z = z) tsplot(sim.cir.a) # simulate from Gallant-Tauchen two factor interest rate diffusion # Set values for the parameters av = .18; avv = -av; avs = -.0088; as = .019; ass = -.0035 b1v = .69; b1vv = 0; b1vs = -.063; b2v = 0; b2s = .038; b2ss = -.017 # Generate the random normals for this simulation n.sim = 10000; n.burn = 10000; ndt = 25 set.seed(0) ird.z = rnorm(2*ndt*(n.sim + n.burn)) # IRD.gensim assumes rho is packed as follows. Note that av is # missing from the list, as the normalization av = -avv is assumed. rho = c(avv, avs, as, ass, b1v, b1vv, b1vs, b2v, b2s, b2ss) ird.sim = IRD.gensim(rho = rho, n.sim = n.sim, n.burn = n.burn, n.var = 2, aux = IRD.aux(ndt = ndt, t.per.sim = 1/52, z = ird.z, X0 = c(1, 5), returnc = c(T, T))) # Plot v(t) and s(t) tsplot(t(matrix(ird.sim, nrow = 2)), lty=c(1,1)) # # illustrate pcode gensim functions # # first simulate OU process using euler1d.pcode.gensim # Generate random normals n.sim = 1000; n.burn = 500; ndt = 25 set.seed(1) z = rnorm(ndt*(n.sim + n.burn)) # set auxiliary parameters ou.names = c("kappa", "theta", "sigma") ou.eu.aux1 = euler.pcode.aux(ndt=25, t.per.sim=1/52, X0 = 0.1, z = z, drift.expr = expression(kappa*(theta - X)), diffuse.expr = expression(sigma), rho.names = ou.names) # test drift and diffusion expressions rho.test = c(0.4, 0.08, 0.1) euler.pcode.test(rho.test,X=0.08,aux=ou.eu.aux1) # simulate using euler1d.pcode.gensim sim.ou.eu1 = euler1d.pcode.gensim(rho=c(0.4, 0.08, 0.1), n.sim = n.sim, n.burn = n.burn, aux = ou.eu.aux1) tsplot(sim.ou.eu1) # truncate lower bound at zero ou.eu.aux1$lbound = 0 sim.ou.eu1 = euler1d.pcode.gensim(rho=c(0.4, 0.08, 0.1), n.sim = n.sim, n.burn = n.burn, aux = ou.eu.aux1) tsplot(sim.ou.eu1) # simulate using euler.pcode.gensim ou.eu.aux = euler.pcode.aux(ndt=25, t.per.sim=1/52, N=1, M=1, X0 = 0.1, z = z, lbound = 0, drift.expr = expression(kappa*(theta - X)), diffuse.expr = expression(sigma), rho.names = ou.names) sim.ou.eu = euler.pcode.gensim(rho = c(kappa, theta, sigma), n.sim = n.sim, n.burn = n.burn, aux = ou.eu.aux1) tsplot(sim.ou.eu) # compare weak 2 with Euler sim.ou.wk = weak2.pcode.gensim(rho = c(kappa, theta, sigma), n.sim = n.sim, n.burn = n.burn, n.var = 1, aux = weak2.pcode.aux(ndt = ndt, t.per.sim = 1/52, N = 1, M = 1, X0 = .1, z = z, drift.expr = expression(kappa*(theta - X)), diffuse.expr = expression(sigma), rho.names = c("kappa", "theta", "sigma"))) # Plot the simulation tsplot(sim.ou.wk) # Plot the absolute differences with the exact solution. tsplot(abs(sim.ou.wk - sim.ou)) lines(abs(sim.ou.eu - sim.ou), col = 2) # multivariate example # Here are the random normals that IRD.gensim used. n.sim = 1000; n.burn = 100; ndt = 25; t.per.sim = 1/52 nsteps = ndt*(n.sim + n.burn) set.seed(0) ird.z = rnorm(2*nsteps) # Make that a matrix with 2 rows, and paste on 4*p + 2 more rows of # new normals. set.seed(5) p = ceiling(.05/(t.per.sim/ndt)) ird.str.z = rbind(matrix(ird.z, nrow = 2), matrix(rnorm((4*p + 2)*nsteps), nrow = 4*p + 2)) # Set values for the parameters av = .18; avv = -av; avs = -.0088; as = .019; ass = -.0035 b1v = .69; b1vv = 0; b1vs = -.063; b2v = 0; b2s = .038; b2ss = -.017 # Pack all parameters into a single vector. Note we include av and avv. rho = c(av, avv, avs, as, ass, b1v, b1vv, b1vs, b2v, b2s, b2ss) ird.sim.str = strong1.pcode.gensim(rho = rho, n.sim = n.sim, n.burn = n.burn, n.var = 2, aux = strong1.pcode.aux(ndt = ndt, t.per.sim = t.per.sim, N = 2, M = 2, X0 = c(1, 5), p = p, z = ird.str.z, drift.expr = expression(av + avv*X[1] + avs*X[2], as + ass*X[2]), diffuse.expr = expression(b1v + b1vv*X[1] + b1vs*X[2], b2v, 0, (b2s + b2ss*X[2])*exp(X[1])), rho.names = c("av","avv","avs", "as", "ass", "b1v", "b1vv", "b1vs", "b2v", "b2s", "b2ss"))) tsplot(t(matrix(ird.sim.str, nrow = 2)),lty=c(1,1)) # Get random normal string for a refinement of the same sample path # as that defined by ird.z. set.seed(3) z.ird.long = refineZ(ird.z, zdim = 2, k = 7) ndt.long = ndt * 128 # 2^7 rho = c(avv, avs, as, ass, b1v, b1vv, b1vs, b2v, b2s, b2ss) ird.sim.ex = IRD.gensim(rho = rho, n.sim = n.sim, n.burn = n.burn, n.var = 2, aux = IRD.aux(ndt = ndt.long, t.per.sim = 1/52, z = z.ird.long, X0 = c(1, 5), returnc = c(T, T))) # Plot the differences with exact between Euler's method and the # strong order 1 scheme. par(mfrow=c(2,1)) tsplot(t(matrix(abs(ird.sim.ex - ird.sim), nrow = 2))) tsplot(t(matrix(abs(ird.sim.ex - ird.sim.str), nrow = 2))) ############################################################################ # simulate from IRD model # ref: Chan, Karolyi, Longstaff and Sanders (1992) - CKLS # An Empirical Comparison of Alternative Models of the Short-Term Interest Rate # Journal of Finance # data: CRSP monthly 30 day t-bill rate ############################################################################ # use CRSP 30 day T-bill data in FinMetrics timeSeries rf.30day # to match data used in CKLS # smpl = (positions(rf.30day) < timeDate("Jan 1990",in.format="%m %Y") # & positions(rf.30day) >= timeDate("June 1964",in.format="%m %Y")) smpl = timeEvent("6/1/1964","1/1/1990") # compute annualized short-rate data shortRate.ts = log(1+rf.30day[smpl])*12 par(mfrow=c(2,1)) plot(shortRate.ts,reference.grid=F,main="") tmp=acf(shortRate.ts) par(mfrow=c(1,1)) summaryStats(shortRate.ts) # GCIR model: dr(t) = (alpha + beta*r(t))dt + sigma*r(t)^gamma *dW(t) # drift = alpha + beta*r(t) # diffusion = sigma*r(t)^gamma # CKLS estimates for 30-day T-bill (annualized rate) sampled monthly # are given in table III pg. 1218 # use euler1d.pcode.gensim to create simulation based on Euler's method n.sim.ckls = 306; n.burn.ckls = 0; ndt.ckls = 25 set.seed(123) z.ckls = rnorm(ndt*(n.sim.ckls + n.burn.ckls)) gcir.names = c("alpha","beta","sigma","gamma") gcir.aux=euler.pcode.aux(drift.expr = expression(alpha + beta*X), diffuse.expr = expression(sigma*X^gamma), rho.names = gcir.names, t.per.sim = 1/12, ndt = 25, z = z.ckls, X0 = 0.067) # gcir model rho.gcir = c(0.0408,-0.5921,sqrt(1.6704),1.4999) gcir.sim = euler1d.pcode.gensim(rho.gcir,n.sim = n.sim.ckls, n.burn = n.burn.ckls, aux = gcir.aux) # merton model model # dr(t) = alpha*dt + sigma*dW(t) rho.merton = c(0.0055,0.0,sqrt(0.0004),0.0) merton.sim = euler1d.pcode.gensim(rho.merton,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # vasicek model model # dr(t) = (alpha+beta*r(t))*dt + sigma*dW(t) rho.vasicek = c(0.0154,-0.1779,sqrt(0.0004),0.0) vasicek.sim = euler1d.pcode.gensim(rho.vasicek,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # cir.sr model model # dr(t) = (alpha+beta*r(t))*dt + sigma*r(t)^0.5*dW(t) rho.cir.sr = c(0.0189,-0.2339,sqrt(0.0073),0.5) cir.sr.sim = euler1d.pcode.gensim(rho.cir.sr,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # dothan model model # dr(t) = sigma*dW(t) rho.dothan = c(0.0,0.0,sqrt(0.1172),1) dothan.sim = euler1d.pcode.gensim(rho.dothan,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # gbm model model # dr(t) = (beta*r(t))*dt + sigma*r(t)*dW(t) rho.gbm = c(0.0,0.1101,sqrt(0.1185),1) gbm.sim = euler1d.pcode.gensim(rho.gbm,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # bs model model # dr(t) = (alpha+beta*r(t))*dt + sigma*r(t)*dW(t) rho.bs = c(0.0242,-0.3142,sqrt(0.1185),1) bs.sim = euler1d.pcode.gensim(rho.bs,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # cir.vr model model # dr(t) = sigma*r(t)^1.5*dW(t) rho.cir.vr = c(0.0,0.0,sqrt(1.5778),1.5) cir.vr.sim = euler1d.pcode.gensim(rho.cir.vr,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) # cev model model # dr(t) = (beta*r(t))*dt + sigma*r(t)^gamma*dW(t) rho.cev = c(0.0,0.1026,sqrt(0.5207),1.2795) cev.sim = euler1d.pcode.gensim(rho.cev,n.sim = n.sim.ckls, n.burn = n.burn.ckls,aux=gcir.aux) par(mfrow=c(3,3)) tsplot(gcir.sim,main="GCIR Model") tsplot(merton.sim,main="Merton Model") tsplot(vasicek.sim,main="Vasicek Model") tsplot(cir.sr.sim,main="CIR SR Model") tsplot(dothan.sim,main="Dothan Model") tsplot(gbm.sim,main="GBM Model") tsplot(bs.sim,main="Brennan-Schwartz Model") tsplot(cir.vr.sim,main="CIR VR Model") tsplot(cev.sim,main="CEV Model") par(mfrow=c(1,1)) ############################################################################ # simulate from 2 factor IRD model # ref Andersen and Lund (1997) # "Estimating continuous-time stochastic volatility models of the short-term # interest rate" # Journal of Econometrics # data: weekly observations (each Wed) of the annualized yield on U.S. T-bills # with 3 months to maturity, constructed from daily series available # from the Fed and corresponds to the 3-month T-bill rates in the weekly # H.15 release. Data are available from 1/61954 through 4/19/1995 # for a total of 2155 observations ############################################################################ # 3-month t-bill data from Gallant and Tauchen are close to data used # by Andersen and Lund colIds(tbill.dat) # [1] "tb3mo" "tb12mo" "tb10yr" plot(tbill.dat[,"tb3mo"],reference.grid=F) summaryStats(tbill.dat[,"tb3mo"]) # general model # dr(t) = k1*(mu - r(t))dt + sig(t)*r(t)^gamma*dW1(t) # dlog(sig2(t)) = k2*(alpha - log(sig2(t))dt + xi*dW2(t) # where gamma > 0 # reparameterized model # dr(t) = k1*(mu - r(t))dt + exp(v(t)/2)*r(t)^gamma*dW1(t) # dv(t) = k2*(alpha - v(t))dt + xi*dW2(t) # k1=0.163, alpha=-0.282, k2=1.04, xi=1.27, mu=5.95, gamma=0.544 # ref: Table 6 rho.AL2.names = c("k1","alpha","k2","xi","mu","gamma") # specify z externally ndt = 25; t.per.sim = 1/52; n.var=2; n.sim=2000; n.burn=100 set.seed(0) z.euler = rnorm(ndt*n.var*(n.sim + n.burn)) ird.AL2.aux = euler.pcode.aux( N = n.var, M = n.var, t.per.sim = t.per.sim, ndt = ndt, X0 = c(6,-0.3), z = z.euler, returnc = c(T,T), rho.names=rho.AL2.names, drift.expr=expression( k1*(mu - X[1]), k2*(alpha-X[2])), diffuse.expr=expression( (exp(X[2]/2))*X[1]^gamma, 0.0, 0.0, xi)) # test drift and diffusion expressions rho.AL2 = c(0.163,-0.282,1.04,1.27,5.95,0.544) euler.pcode.test(rho.AL2, N=2, M=2, t=0, X=c(6,-0.3), aux=ird.AL2.aux) # create simulation and return simulated r(t) and v(t) ird.AL2.sim = euler.pcode.gensim(rho=rho.AL2, n.sim = n.sim, n.var = n.var, n.burn = n.burn, aux = ird.AL2.aux) class(ird.AL2.sim) length(ird.AL2.sim) # reshape data ird.AL2.sim = matrix(ird.AL2.sim,2000,2,byrow=T) par(mfrow=c(2,1)) tsplot(ird.AL2.sim[,1], main="Simulated short rate", ylab="% per year") tsplot(ird.AL2.sim[,2], main="Simulated log-volatility") abline(h=0) par(mfrow=c(1,1)) # simulate using Strong 1 # need z to have 2*ndt*M*(p + 1)*(n.burn+n.sim) elements p.strong = ceiling(0.05/(t.per.sim/ndt)) p.strong set.seed(0) z.strong = rnorm(2*ndt*n.var*(p.strong + 1)*(n.sim + n.burn)) ird.AL2.aux = strong1.pcode.aux( N = n.var, M = n.var, t.per.sim = t.per.sim, ndt = ndt, p = p.strong, z = z.strong, X0 = c(6,-0.3), returnc=c(T,T), rho.names = rho.AL2.names, drift.expr = expression( k1*(mu - X[1]), k2*(alpha-X[2])), diffuse.expr = expression( (exp(X[2]/2))*X[1]^gamma, 0.0, 0.0, xi)) ird.AL2.simS = strong1.pcode.gensim(rho=rho.AL2, n.sim=2000, n.var=2, n.burn=100, aux=ird.AL2.aux) ird.AL2.simS = matrix(ird.AL2.simS,2000,2,byrow=T) par(mfrow=c(2,1)) tsplot(ird.AL2.simS[,1], main="Simulated short rate", ylab="% per year") tsplot(ird.AL2.simS[,2], main="Simulated log-volatility") abline(h=0) par(mfrow=c(1,1)) # simulate using Weak 2 # need z to have (n.burn+n.sim)*ndt*M elements # need u to have (n.burn+n.sim)*ndt*M*(M-1)/2 set.seed(0) z.weak = rnorm(ndt*n.var*(n.sim + n.burn)) u.weak = runif((n.burn+n.sim)*ndt*n.var*(n.var-1)/2) ird.AL2.aux = weak2.pcode.aux( N = n.var, M = n.var, t.per.sim = t.per.sim, ndt = ndt, z = z.weak, u = u.weak, X0 = c(6,-0.3), returnc=c(T,T), rho.names = rho.AL2.names, drift.expr = expression( k1*(mu - X[1]), k2*(alpha-X[2])), diffuse.expr = expression( (exp(X[2]/2))*X[1]^gamma, 0.0, 0.0, xi)) ird.AL2.simW = weak2.pcode.gensim(rho=rho.AL2, n.sim=2000, n.var=2, n.burn=100, aux=ird.AL2.aux) ird.AL2.simW = matrix(ird.AL2.simW,2000,2,byrow=T) par(mfrow=c(2,1)) tsplot(ird.AL2.simW[,1], main="Simulated short rate", ylab="% per year") tsplot(ird.AL2.simW[,2], main="Simulated log-volatility") abline(h=0) par(mfrow=c(1,1)) # put strong 1 and weak 2 plots on 4 panel graph par(mfrow=c(2,2)) tsplot(ird.AL2.simS[,1], main="Simulated short rate from Strong 1", ylab="% per year") tsplot(ird.AL2.simS[,2], main="Simulated log-volatility from Strong 1") abline(h=0) tsplot(ird.AL2.simW[,1], main="Simulated short rate from Weak 2", ylab="% per year") tsplot(ird.AL2.simW[,2], main="Simulated log-volatility from Weak 2") abline(h=0) par(mfrow=c(1,1)) # one factor model # dr(t) = k1*(mu - r(t))dt + sigma*r(t)^gamma * dW(t) # k1=0.082, mu=6.279, sigma=0.670, gamma=0.676 # ref: Table 7 # simulate using Euler # set t.per.sim=1/52 for weekly data # Andersen and Lund set ndt=25 rho.AL1 = c(0.082,6.279,0.670,0.676) rho.AL1.names = c("k1","mu","sigma","gamma") ird.AL1.aux = euler.pcode.aux(drift.expr=expression(k1*(mu-X)), diffuse.expr=expression(sigma*X^gamma), rho.names = rho.AL1.names, t.per.sim = 1/52, ndt = 25, X0 = 6, lbound=0,ubound=100) # check drift and diffusion terms euler.pcode.test(rho.AL1,X=6,aux=ird.AL1.aux) # create simulation ird.AL1.sim = euler1d.pcode.gensim(rho.AL1,n.sim=2000,n.var=1,n.burn=1000, aux=ird.AL1.aux) tsplot(ird.AL1.sim) summaryStats(ird.AL1.sim) par(mfrow=c(2,1)) tmp = acf(ird.AL1.sim) tmp = acf(ird.AL1.sim^2) par(mfrow=c(1,1)) ############################################################################ # simulate from 2 factor IRD model # ref Gallant and Tauchen 2001 # EMM User's Guide Version 1.6 section 7 # ref Gallant and Tauchen 1998 # "Reprojecting partially observed systems with applications to interest # rate diffusions" # JASA # data: 1735 weekly obvs on 3 month T-bill over period 1/5/62 - 3/31/1995 # in timeSeries tbill.dat ############################################################################ # actual interest rate data colIds(tbill.dat) # [1] "tb3mo" "tb12mo" "tb10yr" plot(tbill.dat[,"tb3mo"],reference.grid=F) summaryStats(tbill.dat[,"tb3mo"]) # simple model # dv(t) = 0 # ds(t) = (as+ass*s(t))dt + b2s*exp(v0)dW2(t) # v0 = 1 # as=12.586, ass=-1.264, b2s=0.1176 # calibrated to weekly 3 month t-bill data - annualized rate * 100 # simulatue using Euler's method # note: long-run mean is about 6% # set t.per.sim=1/52 for weekly data # parameters from EMM using short simulation # rho.ird1 = c(12.586,-1.264,0.1176,1) # parameters from EMM using long simulation rho.ird1 = c(6.496,-0.648,-0.0997) rho1.names = c("as","ass","b2s") ird1.aux = euler.pcode.aux(drift.expr=expression(as + ass*X), diffuse.expr=expression(b2s*exp(1)), rho.names = rho1.names, t.per.sim = 1/52, ndt = 25, X0 = 6, lbound=-100,ubound=100) # check drift and diffusion terms euler.pcode.test(rho.ird1,X=6,aux=ird1.aux) # create simulation ird1.sim = euler1d.pcode.gensim(rho.ird1,n.sim=2000,n.var=1,n.burn=500, aux=ird1.aux) tsplot(ird1.sim) # note: mean value looks too high, and SD value is too small summaryStats(ird1.sim) par(mfrow=c(2,1)) tmp = acf(ird1.sim) tmp = acf(ird1.sim^2) par(mfrow=c(1,1)) # general model # dv(t) = (av+avv*v(t)+avs*s(t))dt + (b1v+b1vs*s(t))dW1(t) # ds(t) = (as+ass*s(t))dt + (b2s+b2ss*s(t))*exp(v(t))dW2(t) # av = -avv is normalization used for identification # avv=-0.18375, avs=-0.008245, as=0.01752, ass=-0.003573 # b1v=0.7212, b1vv=0,b1vs=-0.0681, b2v=0.0,b2s=0.03849, b2ss=-0.01692 # calibrated to weekly 3 month t-bill data - annualized rate*100 # use Euler.pcode.gensim rho.ird2 = c(-0.18375,-0.008245,0.01752,-0.003573,0.7212, 0.0,-0.0681,0.0,0.03849,-0.01692) rho2.names = c("avv", "avs", "as", "ass", "b1v", "b1vv", "b1vs", "b2v", "b2s", "b2ss") ird2.aux=euler.pcode.aux( N=2, M=2, ndt = 100, seed=0, lbound=-100, ubound=100, returnc=c(F,T), rho.names=rho2.names, drift.expr=expression( -avv+avv*X[1]+avs*X[2], as+ass*X[2]), diffuse.expr=expression( b1v+b1vv*X[1]+b1vs*X[2], b2v, 0.0, (b2s+b2ss*X[2])*exp(X[1]))) # test drift and diffusion expressions euler.pcode.test(rho.ird2,N=2,M=2,X=c(1,6),aux=ird2.aux) # create simulation and return only simulates s(t) ird2.sim = euler.pcode.gensim(rho=rho.ird2, n.sim=2000, n.var=1, n.burn=100, aux=ird2.aux) tsplot(ird2.sim) summaryStats(ird2.sim) # check simulation with ird.gensim # create simulation and return both s(t) and v(t) ird2.aux$returnc=c(T,T) ird22.sim = euler.pcode.gensim(rho=rho.ird2, n.sim=2000, n.var=2, n.burn=100, aux=ird2.aux) # note - both simulated values are put in one long vector tmp = matrix(ird22.sim,2000,2,byrow=T) par(mfrow=c(2,1)) tsplot(tmp[,1]) tsplot(tmp[,2]) par(mfrow=c(1,1)) summaryStats(tmp) ############################################################################ # simulate from 2 factor continuous time stochastic volatility model # ref Gallant and Tauchen 2001 # EMM User's Guide Version 1.6 section 8 # data: 3778 daily log returns on microsoft stock, adjusted for stock splits # over the period march 13, 1986, through February 23, 2001 # source: ftp: ############################################################################ # general model # dU1(t) = (a10 + a12*U2(t))dt + exp(b10 + b12*U2(t))dW1(t) # dU2(t) = (a20 + a22*U2(t))dt + b20*dW2(t) # a20=0 and b20=1 for identification # a12=0 to eliminate volatility in mean effects # U1(t) = log price process # U2(t) = volatility factor # a10=0.4351, a22=-0.3026, b10=-0.1136, b12=0.2641 # use Euler.pcode.gensim rho.sv2 = c(0.4351,-0.3026,-0.1136,0.2641) rho.sv2.names = c("a10","a22","b10","b12") sv2.aux=euler.pcode.aux( N=2, M=2, t.per.sim=1/252, ndt = 24, seed=0, lbound=-100, ubound=100, X0=c(0,0), returnc=c(T,T), rho.names=rho.sv2.names, drift.expr=expression( a10, a22*X[2]), diffuse.expr=expression( exp(b10+b12*X[2]), 0.0, 0.0, 1)) # test drift and diffusion expressions euler.pcode.test(rho.sv2,N=2,M=2,X=c(0,1),aux=sv2.aux) # create simulation # note: simulated data include log prices and log volatility sv2.sim = euler.pcode.gensim(rho=rho.sv2, n.sim=3800, n.var=2, n.burn=500, aux=sv2.aux) sv2.sim = matrix(sv2.sim,3800,2,byrow=T) par(mfrow=c(2,1)) tsplot(sv2.sim[,1]) tsplot(sv2.sim[,2]) par(mfrow=c(1,1)) # compute log returns ret.sim = diff(sv2.sim[,1])*100 tsplot(ret.sim) summaryStats(sv2.sim) ############################################################################ # simulate from 3 factor continuous time stochastic volatility model # ref Gallant and Tauchen 2002 # Efficient Method of Moments # Working Paper # ############################################################################