# EMMexamples.ssc script file with EMM examples for documentation # author: Eric Zivot # created: January 18, 2004 # updated: April 28, 2004 # updated: April 1, 2005 # # requires S+FinMetrics 2.0 ################################################################################### # MA(1) model # ref: R. Chumacero (1997), "Finite sample properties of the efficient # method of moments," Studies in Nonlinear Dynamics and Econometrics ################################################################################### # simulate from MA(1) using arima.sim with rho=(0,0.5,1) set.seed(123) ma1.sim = arima.sim(model=list(ma=-0.5),n=250) par(mfrow=c(3,1)) tsplot(ma1.sim) abline(h=0) tmp = acf(ma1.sim) tmp = acf(ma1.sim,type="partial") par(mfrow=c(1,1)) # create gensim functions for MA(1) model MA1.gensim <- function(rho, n.sim, n.var=1, n.burn, aux=MA1.aux()) { # simulate from MA(1) model # y(t) = mu + e(t) - psi*e(t-1), e(t) ~ iid N(0,sigma^2) # rho = (mu,psi,sigma2)' # aux is a list with components # z = vector of N(0,1) values of length (n.sim + n.burn) nz <- n.burn + n.sim if(is.null(z <- aux$z)) stop("aux$z must be supplied") if(length(z) != nz) stop("aux$z must be of length",nz," if provided") ans <- rho[1]+arima.sim(model = list(ma=-rho[2]), innov = z[(n.burn+1):nz]*sqrt(rho[3]), start.innov = z[1:n.burn]*sqrt(rho[3])) as.numeric(ans) } # estimate AR(3) auxiliary model using SNP ar3.fit = SNP(data=ma1.sim,model=SNP.model(ar=3)) class(ar3.fit) summary(ar3.fit) plot(ar3.fit,which.plots=3:6) # simulate from estimated AR(3) model ar3.sim = simulate(ar3.fit) par(mfrow=c(3,1)) tsplot(ar3.sim) abline(h=0) tmp = acf(ar3.sim) tmp = acf(ar3.sim,type="partial") par(mfrow=c(1,1)) # fit MA(1) using EMM set.seed(345) nsim = 10000 nburn = 10 MA1.aux = list(z = rnorm(nsim+nburn)) start.vals = c(0,0.5,1) names(start.vals) = c("mu","theta","sigma2") EMM.MA1.fit = EMM(ar3.fit,coef=start.vals, control=EMM.control(n.sim=nsim,n.burn=nburn), gensim.fn="MA1.gensim",gensim.language="SPLUS", gensim.aux=MA1.aux, save.sim=T) names(EMM.MA1.fit$coef) = c("mu","psi","sigma2") class(EMM.MA1.fit) EMM.MA1.fit summary(EMM.MA1.fit) # estimate MA(1) using conditional mle mle.fit = arima.mle(ma1.sim,model=list(ma=-0.5), xreg=rep(1,length(ma1.sim))) mle.fit mle.fit$sigma2 sqrt(mle.fit$var.coef) # fit MA(1) with AR(2) auxiliary model set.seed(123) ar1.sim = arima.sim(model=list(ar=0.5),n=250) ar2.fit2 = SNP(data=as.matrix(ar1.sim),model=SNP.model(ar=2)) summary(ar2.fit2) EMM.fit2 = EMM(ar2.fit2,coef=start.vals, control=EMM.control(n.sim=nsim,n.burn=nburn), gensim.fn="MA1.gensim",gensim.language="SPLUS", gensim.aux=MA1.aux, save.sim=T) names(EMM.fit2$coef) = c("mu","psi","sigma2") EMM.fit2 # note: smdat component of EMM.fit is the simulated data at the optimum # and not the innovations used to simulate the data # EMM diagnostics # note: names for structural model coefs are not printed # and name for variance component is not printed. summary(EMM.fit) ########################################################################################## # discrete time stochastic volatility model ########################################################################################## # Ex 1 # ref: Gallant and Tauchen (2002) # EMM User's Guide Version 1.6 section 3 # data 1: sp500.dat. 4044 daily cc returns * 1000 on S&P 500 index over # the period 1/04/77 - 12/31/92. This data is used in GT (2002) # data 2: sp500.ts. 4365 daily cc returns on S&P 500 index over the period # 3/14/86 - 6/30/2003. Similar data is analyzed in GT (2002) # start(sp500.ts) end(sp500.ts) summaryStats(sp500.ts) plot(sp500.ts,reference.grid=F) abline(h=0) # SNP models # fit.sp500.11118000 best BIC manual search # fit.sp500.auto best BIC from SNP.auto # fit.sp500.20b14000 preferred model from GT (using older data) # fit.sp500.20b14010 preferred model from GT (using older data) fit.sp500.11118000 # plot standardized residuals, SNP density and ACF of squared std residuals plot(fit.sp500.11118000,ask=F,which.plots=c(2,3,7)) # # simulator function for GT's simple discrete time SV model # # y(t) - mu = c(y(t-1) - mu) + sigma*exp(w(t))*z(t) # w(t) = sigma.w*(zw(t) - sigma.w) # sv.gt.gensim <- function(rho,n.sim,n.var=1,n.burn,aux=sv.gt.aux) { # rho = (mu,c,sigma,sigma.w) # aux is a list with components # z = vector of N(0,1) values of length (n.sim + n.burn) # zw = vector N(0,1) values of length (n.sim + n.burn) nz = n.burn + n.sim sv = rho[3]*exp(-rho[4]*rho[4] + rho[4]*aux$zw)*aux$z y = arima.sim(model = list(ar=rho[2]), innov = sv[(n.burn+1):nz], start.innov = sv[1:n.burn]) y = y + rho[1] as.numeric(y) } # simulate SV model calibrated to match sp500.ts n.sim=4000 n.burn=100 nz = n.sim+n.burn set.seed(123) sv.gt.aux = list(z=rnorm(nz),zw=rnorm(nz)) rho.gt = c(0.0004,0.01,0.01,.5) names(rho.gt) = c("mu","c","sigma","sigma.w") sv.gt.sim = sv.gt.gensim(rho=rho.gt,n.sim=n.sim,n.burn=n.burn,aux=sv.gt.aux) tsplot(sv.gt.sim) abline(h=0) summaryStats(sv.gt.sim) # fit SNP model to simulated data fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) fit.sp500.auto = SNP.auto(as.matrix(sv.gt.sim),n.drop=14, control = SNP.control(xTransform="logistic", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), arMax=1,xPolyMax=1,lagPMax=4) # EMM estimation using SNP 10014000 auxiliary model n.sim = 10000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.gt.aux = list(z=rnorm(nz),zw=rnorm(nz)) rho.gt = c(0.0004,0.01,0.01,.5) names(rho.gt) = c("mu","c","sigma","sigma.w") emm.svgt.fit1 = EMM(fit.sp500.10014000,coef=rho.gt, EMM.control(n.burn=n.burn,n.sim=n.sim), gensim.fn="sv.gt.gensim",gensim.language="SPLUS", gensim.aux=sv.gt.aux) names(emm.svgt.fit1$coef) = names(rho.gt) summary(emm.svgt.fit1) # EMM estimation using SNP 11118000 auxiliary model n.sim = 10000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.gt.aux = list(z=rnorm(nz),zw=rnorm(nz)) rho.gt = c(0.0004,0.01,0.01,.5) names(rho.gt) = c("mu","c","sigma","sigma.w") emm.svgt.fit = EMM(fit.sp500.11118000,coef=rho.gt, EMM.control(n.burn=n.burn,n.sim=n.sim), gensim.fn="sv.gt.gensim",gensim.language="SPLUS", gensim.aux=sv.gt.aux) names(emm.svgt.fit$coef) = names(rho.gt) summary(emm.svgt.fit) # Ex. 2 # Discrete time stochastic volatility model # ref. Andersen and Sorensen (1996) # "GMM Estimation of a Stochastic Volatility Model: A Monte Carlo Study", # Journal of Business and Economic Statistics # ref. Andersen, Chung and Sorensen (1999) # "Efficient Method of Moments Estimation of a Stochastic Volatility # Model: A Monte Carlo Study" # Journal of Econometrics, 91. # # # simulator function for AS SV model # # y(t) = exp(w(t)/2)*z(t) # w(t) = alpha + beta*w(t-1) + sigma.u*u(t) # -1 < beta < 1, sigma.u > 0 # # simulator function # note: need transformation to impose stationarity on simulations # set beta = exp(beta.t)/(1+exp(beta.t)) # beta.t = log(beta/(1-beta)) sv.as.gensim <- function(rho,n.sim,n.var=1,n.burn,aux=sv.as.aux) { # rho = (alpha, beta.t, sigma.u) # aux is a list with components # z = vector of N(0,1) values of length (n.sim + n.burn) # u = vector N(0,1) values of length (n.sim + n.burn) nz = n.burn + n.sim beta = exp(rho[2])/(1+exp(rho[2])) mu = rho[1]/(1-beta) w = mu + arima.sim(model = list(ar=beta), innov = aux$u[(n.burn+1):nz]*rho[3], start.innov = aux$u[1:n.burn]*rho[3]) y = exp(w/2)*aux$z[(n.burn+1):nz] as.numeric(y) } # simulate data from AS Monte Carlo design 1 # alpha=-0.736, beta=0.90, sigma.u=0.363 # beta.t = log(beta/(1-beta)) n.sim = 1000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.as.aux = list(z=rnorm(nz),u=rnorm(nz)) rho.as = c(-0.736,2.197225,0.363) names(rho.as) = c("alpha","beta.t","sigma.u") sv.as.sim = sv.as.gensim(rho=rho.as,n.sim=n.sim, n.burn=n.burn,aux=sv.as.aux) tsplot(sv.as.sim) abline(h=0) summaryStats(sv.as.sim) # estimate GARCH(1,1) auxiliary model to simulated data fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) fit.svassim.01110000 <- SNP(as.matrix(sv.as.sim), model=SNP.model(ar=0,arch=1,garch=1), control = SNP.control(xTransform="logistic", n.start=n.start, fOld=fOld, fNew=fNew), n.drop=4, trace=T) # estimate by EMM n.sim = 10000 n.burn = 100 nz = n.sim+n.burn set.seed(123) sv.as.aux = list(z=rnorm(nz),u=rnorm(nz)) emm.svassim.fit = EMM(fit.svassim.01110000,coef=rho.as, EMM.control(n.burn=n.burn,n.sim=n.sim), gensim.fn="sv.as.gensim",gensim.language="SPLUS", gensim.aux=sv.as.aux) names(emm.svassim.fit$coef) = names(rho.as) # alpha=-0.736, beta=0.90, sigma.u=0.363 summary(emm.svassim.fit) # simulate data from AS Monte Carlo design 2 # calibrated to daily data # alpha=-0.147, beta=0.98, sigma.u=0.166 # beta.t = log(beta/(1-beta)) n.sim = 4000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.as2.aux = list(z=rnorm(nz),u=rnorm(nz)) rho.as2 = c(-0.147,3.89182,0.166) names(rho.as2) = c("alpha","beta.t","sigma.u") sv.as2.sim = sv.as.gensim(rho=rho.as2,n.sim=n.sim, n.burn=n.burn,aux=sv.as2.aux) tsplot(sv.as2.sim) abline(h=0) summaryStats(sv.as2.sim) # estimate GARCH(1,1) auxiliary model to simulated data fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) fit.svas2sim.01110000 <- SNP(as.matrix(sv.as2.sim), model=SNP.model(ar=0,arch=1,garch=1), control = SNP.control(xTransform="logistic", n.start=n.start, fOld=fOld, fNew=fNew), n.drop=4, trace=T) # estimate by EMM n.sim = 40000 n.burn = 100 nz = n.sim+n.burn set.seed(123) sv.as2.aux = list(z=rnorm(nz),u=rnorm(nz)) emm.svas2sim.fit = EMM(fit.svas2sim.01110000,coef=rho.as2, EMM.control(n.burn=n.burn,n.sim=n.sim), gensim.fn="sv.as.gensim",gensim.language="SPLUS", gensim.aux=sv.as2.aux) names(emm.svas2sim.fit$coef) = names(rho.as2) summary(emm.svas2sim.fit) # delta method computation of standard error bt.hat = emm.svas2sim.fit$coef[2] b.hat = exp(bt.hat)/(1+exp(bt.hat)) dg = b.hat*(1 - b.hat) se.b.hat = sqrt(dg%*%emm.svas2sim.fit$var[2,2]%*%dg) se.b.hat # compare with GMM estimation # create data to be passed to sv.moments function sv.pow = cbind(abs(sv.as2.sim),sv.as2.sim^2,abs(sv.as2.sim)^3, sv.as2.sim^4) sv2.pow = sv.pow[-(1:10),] sv2.cm = tslag(sv.as2.sim, 1:10, trim=T) * as.vector(sv.as2.sim[-(1:10)]) sv2.data = cbind(sv2.pow, abs(sv2.cm), sv2.cm^2) colIds(sv2.data) = NULL # iterative GMM estimation of the SV model start.vals = c(0,0.5,0.5) names(start.vals) = c("omega","beta","sigu") sv2.fit.gmm = GMM(start.vals, sv.moments, method="iterative", ts=T, data=sv2.data) # note: results are remarkably close to EMM estimates! summary(sv2.fit.gmm) # estimate SV model using sp500 data n.sim = 40000 n.burn = 100 nz = n.sim+n.burn set.seed(456) sv.as.aux = list(z=rnorm(nz),u=rnorm(nz)) rho.as = c(-0.147,3.89182,0.166) names(rho.as) = c("alpha","beta.t","sigma.u") emm.svas.sp500 = EMM(fit.sp500.11118000,coef=rho.as, EMM.control(n.burn=n.burn,n.sim=n.sim), gensim.fn="sv.as.gensim",gensim.language="SPLUS", gensim.aux=sv.as2.aux) names(emm.svas.sp500$coef) = names(rho.as) summary(emm.svas.sp500) # # Ex 3 General univariate discrete-time SV model # ref. Gallant, Hsieh and Tauchen (1997), "Estimation of stochastic # volatility models with diagnostics," Journal of Econometrics # 81, 159-192 # ?SVOL args(SVOL) # estimate AS SV model using SVOL function set.seed(456) start.vals = c(0,exp(-0.4/2), 0.9, 0.1) svol.fit1 <- SVOL(fit.svassim.01110000, order=c(0,1), param.ini=start.vals, gensim.language="C", control=EMM.control(n.sim=40000, n.burn=100)) ########################################################################### # continuous time interest rate models ########################################################################### # # 1. Interest rate diffusion models # ref: Chan, Karolyi, Longstaff and Sanders (1992) # Journal of Finance # James and Webber (2000). Interest Rate Models # # CKLS model: dr(t) = (alpha + beta*r(t))dt + sigma*r(t)^gamma *dW(t) # drift = alpha + beta*r(t) # diffusion = sigma*r(t)^gamma # re-parametrized drift: kappa(mu - r(t)) => alpha = kappa*mu; beta = -kappa # CKLS parameterization for 30-day T-bill (annualized rate) # Table III page 1218 # alpha = 0.0408, beta = -0.5921, sigma2 = 1.6704, gamma = 1.4999 # long-run mean = alpha/-beta = 0.0689 # use euler1d.pcode.gensim to create simulation based on Euler's method ?euler1d.pcode.gensim args(euler1d.pcode.gensim) # function(rho, n.sim, n.var = 1, n.burn = 0, aux = euler.pcode.aux()) ?euler.pcode.aux args(euler.pcode.aux) # function(drift.expr = expression(0), diffuse.expr = expression(0), # rho.names = character(0), N = 1, M = 1, t.per.sim = 1, ndt = 100, # z = NULL, seed = 0, X0 = rep(0.1, length = N), returnc = rep(T, # length = N), lbound = 0, ubound = 100) rho.names = c("alpha","beta","sigma","gamma") ckls.aux = euler.pcode.aux( drift.expr=expression(alpha + beta*X), diffuse.expr=expression(sigma*X^gamma), rho.names = rho.names, t.per.sim = 1/12, ndt = 25, lbound=0,ubound=10) # auxiliary information based on re-parameterized model rho2.names = c("mu","kappa","sigma","gamma") ckls.aux2 = euler.pcode.aux( drift.expr=expression(kappa*(mu - X)), diffuse.expr = expression(sigma*X^gamma), rho.names = rho2.names, t.per.sim = 1/12, ndt = 25, lbound=0, ubound=10) # test the drift and diffusion expressions using CKLS parameter estimates rho2.ckls = c(0.0408/0.5921, 0.5921, sqrt(1.6704), 1.4999) euler.pcode.test(rho2.ckls, X=0.05, aux=ckls.aux2) # create sample simulation based on CKLS estimates # using CKLS parameterization rho.ckls = c(0.0408, -0.5921, sqrt(1.6704), 1.4999) ckls.aux$seed = 123 ckls.aux$X0 = 0.0408/0.5921 ckls.sim = euler1d.pcode.gensim(rho.ckls,n.sim=300, n.burn=100, aux=ckls.aux) class(ckls.sim) tsplot(ckls.sim) summaryStats(ckls.sim) par(mfrow=c(2,1)) tmp = acf(ckls.sim) tmp = acf(ckls.sim^2) par(mfrow=c(1,1)) # simulate from re-parameterized model using euler1d.pcode.gensim rho2.ckls = c(0.0408/0.5921, 0.5921, sqrt(1.6704), 1.4999) ckls.aux2$seed = 123 ckls.aux2$X0 = 0.0408/0.5921 ckls.sim2 = euler1d.pcode.gensim(rho2.ckls, n.sim=300, n.burn=100, aux=ckls.aux2) # check that simulations are the same tsplot(ckls.sim2,ckls.sim) tsplot(ckls.sim2) summaryStats(ckls.sim2) ckls.resid = residuals(OLS(ckls.sim~tslag(ckls.sim),na.rm=T)) tsplot(ckls.resid) par(mfrow=c(3,1)) tsplot(ckls.sim2) tmp = acf(ckls.sim2) tmp = acf(ckls.resid^2) par(mfrow=c(1,1)) summaryStats(ckls.sim2) # automatically fit SNP model to simulated interest rates fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) # fit SNP model ckls.sim = as.matrix(ckls.sim) fit.snp.ckls = SNP.auto(ckls.sim,n.drop=8, control = SNP.control(xTransform="spline", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), xPolyMax=1,lagPMax=1) fit.snp.ckls summary(fit.snp.ckls) # snp diagnostics plot(fit.snp.ckls) snp.ckls.resid = residuals(fit.snp.ckls,standardized=F) par(mfrow=c(3,1)) tsplot(snp.ckls.resid,main="standardized residuals") abline(h=0) tmp = acf(snp.ckls.resid) tmp = acf(snp.ckls.resid^2) par(mfrow=c(1,1)) # simulate from fitted SNP model snp.ckls.sim = simulate(fit.snp.ckls) par(mfrow=c(3,1)) tsplot(snp.ckls.sim) tmp = acf(snp.ckls.sim) tmp = acf(snp.ckls.resid^2) par(mfrow=c(1,1)) summaryStats(snp.ckls.sim) # EMM estimation of ckls model # note: sample size is too small to get stable estimates n.burn = 100 n.sim = 40000 ndt = 25 set.seed(456) z.ckls = rnorm(ndt*(n.burn + n.sim)) rho.ckls.names = c("alpha","beta","sigma","gamma") ckls.aux = euler.pcode.aux( drift.expr = expression(alpha + beta*X), diffuse.expr = expression(sigma*X^gamma), rho.names = rho.ckls.names, t.per.sim = 1/12, ndt = 25, z = z.ckls, X0 = 0.06, lbound=0, ubound=10) rho.ckls = c(0.04, -0.6, 1.3 , 1) emm.ckls.fit = EMM(fit.snp.ckls, coef=rho.ckls, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler1d.pcode.gensim, gensim.language = "SPLUS", gensim.aux = ckls.aux, save.sim=T) names(emm.ckls.fit$coef) = rho.ckls.names emm.ckls.fit # The high objective value is strange. # It looks like it is important to use random re-starts # since there appears to be multiple local minimums # To see this, re-estimate with a different starting value for gamma rho.ckls = c(0.04, -0.6, 1.3 , 1.5) emm.ckls.fit2 = EMM(fit.snp.ckls,coef=rho.ckls, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler1d.pcode.gensim, gensim.language = "SPLUS", gensim.aux=ckls.aux,save.sim=T) names(emm.ckls.fit2$coef) = rho.ckls.names emm.ckls.fit2 # now do estimation with random re-starts rho.ckls = c(0.04, -0.6, 1.3 , 1) emm.ckls.fit.restart = EMM(fit.snp.ckls,coef=rho.ckls, control = EMM.control(n.sim=n.sim,n.burn=n.burn, n.start = rep(5,3), tweak=c(0.25, 0.5, 1.0)), gensim.fn = euler1d.pcode.gensim, gensim.language = "SPLUS", gensim.aux = ckls.aux, save.sim=T) names(emm.ckls.fit.restart$coef) = rho.ckls.names emm.ckls.fit.restart summary(emm.ckls.fit.restart) # # compare to GMM estimation of CKLS model # # function to compute CKLS moments # note: dt defines the discretization step. For example, # if r(t) is annualized rate sampled monthly, set dt = 1/12 # if r(t) is monthly rate sampled monthly, set dt = 1 ckls.moments <- function(parm,data=NULL,dt=1/12) { # parm = (alpha,beta,sigma,gamma)' # data = [r(t+dt)-r(t),r(t)] # dt = discretization step e.hat = as.vector(data[,1] - (parm[1] + parm[2]*data[,2])*dt) m2 = e.hat*as.vector(data[,2]) m3 = e.hat^2 - dt*parm[3]*parm[3]*(as.vector(data[,2])^(2*parm[4])) m4 = m3*data[,2] cbind(e.hat,m2,m3,m4) } data.ckls.sim = cbind(diff(ckls.sim),tslag(ckls.sim,trim=T)) colIds(data.ckls.sim) = c("RF","RF.diff") # GMM estimation of just identified ckls model start.vals = c(0.06, -0.5, 1, 1) names(start.vals) = c("alpha","beta","sigma","gamma") gmm.ckls.sim = GMM(start.vals,ckls.moments, ts=T, data=data.ckls.sim, dt=1/12) summary(gmm.ckls.sim) # # try estimation with much larger sample size # rho.ckls = c(0.0408,-0.5921,sqrt(1.6704),1.4999) rho.names = c("alpha","beta","sigma","gamma") ckls.aux = euler.pcode.aux(drift.expr=expression(alpha + beta*X), diffuse.expr=expression(sigma*X^gamma), rho.names = rho.names, t.per.sim = 1/12, ndt = 25, seed = 123, lbound=0,ubound=10, X0 = 0.0408/0.5921) ckls.sim2 = euler1d.pcode.gensim(rho.ckls,n.sim=1000,n.var=1,n.burn=100, aux=ckls.aux) tsplot(ckls.sim2) summaryStats(ckls.sim2) # EMM estimation of ckls model # note: sample size is too small to get stable estimates n.burn = 100 n.sim = 10000 set.seed(456) z = rnorm(ndt*(n.burn+n.sim)) rho.ckls = c(0.04, -0.6, 1.3 , 1) rho.ckls.names = c("alpha","beta","sigma","gamma") ckls.aux = euler.pcode.aux(drift.expr=expression(alpha + beta*X), diffuse.expr=expression(sigma*X^gamma), rho.names = rho.ckls.names, t.per.sim = 1/12, ndt = 25, z = z, X0 = 0.06, lbound=0,ubound=10) # fit SNP model ckls.sim2 = as.matrix(ckls.sim2) fit.snp2.ckls = SNP.auto(ckls.sim2,n.drop=8, control = SNP.control(xTransform="spline", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), xPolyMax=1,lagPMax=1) emm.ckls.fit3 = EMM(fit.snp2.ckls,coef=rho.ckls, control=EMM.control(n.sim=n.sim,n.burn=n.burn), gensim.fn=euler1d.pcode.gensim,gensim.language="SPLUS", gensim.aux=ckls.aux,save.sim=T) names(emm.ckls.fit3$coef) = rho.ckls.names emm.ckls.fit3 # # Mike Cliff's monthly T-Bill data which is close to the data used # in CKLS # colIds(ckls.ts) start(ckls.ts) end(ckls.ts) plot(ckls.ts,reference.grid=F) summaryStats(ckls.ts) # fit SNP model fit.snp.ckls.ts = SNP.auto(ckls.ts,n.drop=8, control = SNP.control(xTransform="spline", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), xPolyMax=1,lagPMax=1) # diagnostics plot(fit.snp.ckls.ts, ask=F, which.plots=c(5,7)) plot(simulate(fit.snp.ckls.ts), reference.grid=F) # fit model by EMM n.burn = 100 n.sim = 40000 ndt = 25 set.seed(456) z.ckls = rnorm(ndt*(n.burn + n.sim)) rho.ckls.names = c("alpha","beta","sigma","gamma") ckls.aux = euler.pcode.aux( drift.expr = expression(alpha + beta*X), diffuse.expr = expression(sigma*X^gamma), rho.names = rho.ckls.names, t.per.sim = 1/12, ndt = 25, z = z.ckls, X0 = 0.06, lbound=0, ubound=10) rho.ckls = c(0.04, -0.6, 1.3 , 1.5) emm.ckls.ts.fit = EMM(fit.snp.ckls.ts,coef=rho.ckls, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler1d.pcode.gensim, gensim.language = "SPLUS", gensim.aux=ckls.aux,save.sim=T) names(emm.ckls.ts.fit$coef) = rho.ckls.names emm.ckls.ts.fit # # 2. SV models for interest rates # 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/6/1954 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) tb3mo = tbill.dat[,"tb3mo"] plot(tb3mo,main="Weekly 3 Month U.S. T-Bill Rate", reference.grid=F) summaryStats(tb3mo) par(mfrow=c(2,1)) tmp = acf(tb3mo) tmp = acf(tb3mo,type="partial") par(mfrow=c(1,1)) colIds(tbill.dat) # [1] "tb3mo" "tb12mo" "tb10yr" plot(tbill.dat[,"tb3mo"],reference.grid=F) summaryStats(tbill.dat[,"tb3mo"]) # one factor generalized CIR (gcir) model for weekly interest rate data # dr(t) = k*(mu - r(t))dt + sigma*r(t)^gamma * dW(t) # k1=0.082, mu=6.279, sigma=0.670, gamma=0.676 # ref: Andersen and Lund, Table 7 page 370 # simulate generalized CIR process using Euler's method # set t.per.sim=1/52 for weekly data # Andersen and Lund set ndt=25 rho.gcir.names = c("k","mu","sigma","gamma") gcir.aux = euler.pcode.aux( drift.expr = expression(k*(mu-X)), diffuse.expr = expression(sigma*X^gamma), rho.names = rho.gcir.names, t.per.sim = 1/52, ndt = 25, lbound=0,ubound=100) # check drift and diffusion terms euler.pcode.test(rho.gcir,X=6,aux=gcir.aux) # simulate 2000 values rho.gcir = c(0.082,6.279,0.670,0.676) gcir.aux$seed = 123 gcir.aux$X0 = 6 gcir.sim = euler1d.pcode.gensim(rho.gcir, n.sim=2000, n.var=1, n.burn=1000, aux=gcir.aux) tsplot(gcir.sim) summaryStats(gcir.sim) # fit SNP model to simulated GCIR model for weekly interest rates # automatically fit SNP model to simulated interest rates fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) # fit SNP model gcir.sim = as.matrix(gcir.sim) fit.snp.gcir = SNP.auto(gcir.sim,n.drop=8, control = SNP.control(xTransform="spline", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), xPolyMax=1,lagPMax=1) fit.snp.gcir summary(fit.snp.gcir) # snp diagnostics plot(fit.snp.gcir) snp.gcir.resid = residuals(fit.snp.gcir,standardized=T) par(mfrow=c(3,1)) tsplot(snp.gcir.resid,main="standardized residuals") abline(h=c(2,-2)) tmp = acf(snp.gcir.resid) tmp = acf(snp.gcir.resid^2) par(mfrow=c(1,1)) # simulate from fitted SNP model snp.gcir.sim = simulate(fit.snp.gcir) tsplot(snp.gcir.sim,main="simulation from SNP 11111010 model") # EMM estimation of gcir model n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.gcir = rnorm(ndt*(n.burn+n.sim)) rho.gcir = c(0.082,6.279,0.670,0.676) rho.gcir.names = c("k","mu","sigma","gamma") gcir.aux = euler.pcode.aux( drift.expr = expression(k*(mu-X)), diffuse.expr = expression(sigma*X^gamma), rho.names = rho.gcir.names, t.per.sim = 1/52, ndt = 25, z = z.gcir, X0 = 6, lbound=0, ubound=100) emm.gcir.fit = EMM(fit.snp.gcir,coef=rho.gcir, control=EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = gcir.aux, save.sim=T) names(emm.gcir.fit$coef) = rho.gcir.names emm.gcir.fit summary(emm.fit) # # fit one-factor model to weekly observations on US T-Bills # # Use the following SNP models for the fits: # 11118000, 10818000, 51118000, 10414010 # in the time series objects # fit.tb3mo.11118000, fit.tb3mo.10818000, # fit.tb3mo.auto, fit.tb3mo.10414010 # 11118000 model # estimate with nsim = 50000 n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.gcir= rnorm(ndt*(n.burn + n.sim)) gcir.aux$z = z.gcir rho.gcir = c(0.082,6.279,0.670,0.676) emm.gcir.11118000.fit = EMM(fit.tb3mo.11118000,coef=rho.gcir, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = gcir.aux, save.sim=T) names(emm.gcir.11118000.fit$coef) = rho.gcir.names emm.gcir.11118000.fit #10818000 model emm.gcir.10818000.fit = EMM(fit.tb3mo.10818000,coef=rho.gcir, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = gcir.aux, save.sim=T) names(emm.gcir.10818000.fit$coef) = rho.gcir.names emm.gcir.10818000.fit #51118000 model emm.gcir.51118000.fit = EMM(fit.tb3mo.auto,coef=rho.gcir, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = gcir.aux, save.sim=T) names(emm.gcir.51118000.fit$coef) = rho.gcir.names emm.gcir.51118000.fit #10414010 model emm.gcir.10414010.fit = EMM(fit.tb3mo.10414010,coef=rho.gcir, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = gcir.aux, save.sim=T) names(emm.gcir.10414010.fit$coef) = rho.gcir.names emm.gcir.10414010.fit # compare diagnostics for the 4 fits par(mfrow=c(2,2)) plot(emm.gcir.11118000.fit) title(sub="SNP 11118000 Model") plot(emm.gcir.10818000.fit) title(sub="SNP 10818000 Model") plot(emm.gcir.51118000.fit) title(sub="SNP 51118000 Model") plot(emm.gcir.10414010.fit) title(sub="SNP 10414010 Model") par(mfrow=c(1,1)) # general 2 factor 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 # trick: shift v(t) up by 5 and set lower bound at zero rho.tf = c(0.163,5.95,0.544,1.04,5-0.282,1.27) rho.tf.names = c("k1","mu","gamma","k2","alpha","xi") tf.aux=euler.pcode.aux( rho.names=rho.tf.names, drift.expr=expression( k1*(mu - X[1]), k2*(alpha-X[2])), diffuse.expr=expression( (exp((X[2]-5)*0.5))*X[1]^gamma, 0.0, 0.0, xi), N=2, M=2, t.per.sim=1/52, ndt = 25, lbound=0,ubound=100, returnc=c(T,T)) # test drift and diffusion expressions euler.pcode.test(rho.tf, N=2, M=2, X=c(6,5-0.3), aux=tf.aux) # create simulation tf.aux$seed = 678 tf.aux$X0 = c(6,5-0.28) tf.sim = euler.pcode.gensim(rho=rho.tf, n.sim=2000, n.var=2, n.burn=1000, aux=tf.aux) tmp = matrix(tf.sim,2000,2,byrow=T) tfird.sim = tmp[,1, drop=F] logvol.sim = tmp[,2]-5 par(mfrow=c(2,1)) tsplot(tmp[,1],main="interest rates: r(t)") tsplot(tmp[,2]-5,main="log volatility: v(t) - 5") abline(h=0) par(mfrow=c(1,1)) summaryStats(tmp[,1]) summaryStats(tmp[,2]-5) # automatically fit SNP model to simulated interest rates fOld <- c(0,1e-5, 0, -1e-5, 1e-5, 1e-5, 1e-4, 0, -1e-4, 1e-4, 1e-4, 1e-3, 0, -1e-3, 1e-3, 1e-3, 1e-2, 0, -1e-2, 1e-2, 1e-2, 1e-1, 0, -1e-1, 1e-1, 1e-1, 1e0, 0, -1e0, 1e0, 1e0) fNew <- c(0,0, 1e-5, 1e-5, -1e-5, 1e-5, 0, 1e-4, 1e-4, -1e-4, 1e-4, 0, 1e-3, 1e-3, -1e-3, 1e-3, 0, 1e-2, 1e-2, -1e-2, 1e-2, 0, 1e-1, 1e-1, -1e-1, 1e-1, 0, 1e0, 1e0, -1e0, 1e0) n.start <- c(0,rep(25,30)) # fit SNP model fit.snp.tfird = SNP.auto(tfird.sim,n.drop=8, control = SNP.control(xTransform="spline", n.start = n.start, seed = 011667, fOld = fOld, fNew = fNew), lagPMax = 1) fit.snp.tfird summary(fit.snp.tfird) # snp diagnostics plot(fit.snp.tfird) snp.tfird.resid = residuals(fit.snp.tfird,standardized=T) par(mfrow=c(3,1)) tsplot(snp.tfird.resid,main="standardized residuals") abline(h=c(2,-2)) tmp = acf(snp.tfird.resid) tmp = acf(snp.tfird.resid^2) par(mfrow=c(1,1)) # simulate from fitted SNP model set.seed(456) snp.tfird.sim = simulate(fit.snp.tfird) tsplot(snp.tfird.sim,main="simulation from SNP 11114010 model") # plot estimated conditional volatility and actual volatility snp.sigma.t = sigma.t.SNP(fit.snp.tfird) par(mfrow=c(2,1)) tsplot(tfird.sim) tsplot(snp.sigma.t, exp(logvol.sim/2)) par(mfrow=c(1,1)) # EMM estimation of tfird model n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) rho.tfird = c(0.163,5.95,0.544,1.04,5-0.282,1.27) rho.tfird.names = c("k1","mu","gamma","k2","alpha","xi") tfird.aux=euler.pcode.aux( rho.names=rho.tfird.names, drift.expr=expression( k1*(mu - X[1]), k2*(alpha-X[2])), diffuse.expr=expression( (exp((X[2]-5)*0.5))*X[1]^gamma, 0.0, 0.0, xi), N=2, M=2, t.per.sim=1/52, ndt = 25, z = z.tfird, X0 = c(6,5-0.28), lbound=0,ubound=100, returnc=c(T,F)) emm.tfird.fit = EMM(fit.snp.tfird,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.tfird.fit$coef) = rho.tfird.names emm.tfird.fit summary(emm.tfird.fit) # # fit two factor model to weekly interest rate data # Use the following SNP models for the fits: # 11118000, 10818000, 51118000, 10414010 # in the time series objects # fit.tb3mo.11118000, fit.tb3mo.10818000, # fit.tb3mo.auto, fit.tb3mo.10414010 # 11118000 model # estimate with nsim = 50000 n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird rho.tfird = c(0.163,5.95,0.544,1.04,5-0.282,1.27) emm.11118000.fit = EMM(fit.tb3mo.11118000,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.11118000.fit$coef) = rho.tfird.names emm.11118000.fit plot(emm.11118000.fit) # estimate with nsim = 75000 n.burn = 1000 n.sim = 75000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird emm.11118000.fit2 = EMM(fit.tb3mo.11118000,coef=emm.11118000.fit$coef, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.11118000.fit2$coef) = rho.tfird.names emm.11118000.fit2 plot(emm.11118000.fit2) # estimate with nsim = 75000 and nburn = 5000 n.burn = 5000 n.sim = 75000 ndt = 24 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird tfird.aux$ndt = ndt tfird.aux$X0 = c(6, 5-0.3) emm.11118000.fit3 = EMM(fit.tb3mo.11118000,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.11118000.fit3$coef) = rho.tfird.names emm.11118000.fit3 plot(emm.11118000.fit3) # 10818000 model # estimate with nsim = 50000 n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird rho.tfird = c(0.163,5.95,0.544,1.04,5-0.282,1.27) emm.10818000.fit = EMM(fit.tb3mo.10818000,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.10818000.fit$coef) = rho.tfird.names # note: std errors are large and the 95% confidence intervals look # weird emm.10818000.fit plot(emm.10818000.fit) # estimate with nsim = 75000 n.burn = 1000 n.sim = 75000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird emm.10818000.fit2 = EMM(fit.tb3mo.10818000,coef=emm.10818000.fit$coef, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.10818000.fit2$coef) = rho.tfird.names emm.10818000.fit2 plot(emm.10818000.fit2) # try estimation using starting values from 11118000 model # results converge to essentially the same point n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird emm.10818000.fit3 = EMM(fit.tb3mo.10818000,coef=emm.11118000.fit$coef, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.10818000.fit3$coef) = rho.tfird.names emm.10818000.fit3 # 51118000 model # estimate with nsim = 50000 n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird rho.tfird = c(0.163,5.95,0.544,1.04,5-0.282,1.27) emm.51118000.fit = EMM(fit.tb3mo.auto,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.51118000.fit$coef) = rho.tfird.names # note: std errors are large and the 95% confidence intervals look # weird emm.51118000.fit plot(emm.51118000.fit) # estimate with nsim = 75000 n.burn = 1000 n.sim = 75000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird emm.51118000.fit2 = EMM(fit.tb3mo.auto,coef=emm.51118000.fit$coef, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.51118000.fit2$coef) = rho.tfird.names emm.51118000.fit2 plot(emm.51118000.fit2) # 10414010 model # estimate with nsim = 50000 n.burn = 1000 n.sim = 50000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird rho.tfird = c(0.163,5.95,0.544,1.04,5-0.282,1.27) emm.10414010.fit = EMM(fit.tb3mo.10414010,coef=rho.tfird, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.10414010.fit$coef) = rho.tfird.names # note: std errors are large and the 95% confidence intervals look # weird emm.10414010.fit plot(emm.10414010.fit) # estimate with nsim = 75000 n.burn = 1000 n.sim = 75000 ndt = 25 set.seed(456) z.tfird = rnorm(2*ndt*(n.burn+n.sim)) tfird.aux$z = z.tfird emm.10414010.fit2 = EMM(fit.tb3mo.10414010,coef=emm.10414010.fit$coef, control = EMM.control(n.sim=n.sim, n.burn=n.burn), gensim.fn = euler.pcode.gensim, gensim.language = "SPLUS", gensim.aux = tfird.aux, save.sim=T) names(emm.10414010.fit2$coef) = rho.tfird.names emm.10414010.fit2 plot(emm.10414010.fit2) # compare diagnostics for the 4 fits par(mfrow=c(2,2)) plot(emm.11118000.fit2) title(sub="SNP 11118000 Model") plot(emm.10818000.fit2) title(sub="SNP 10818000 Model") plot(emm.51118000.fit2) title(sub="SNP 51118000 Model") plot(emm.10414010.fit2) title(sub="SNP 10414010 Model") par(mfrow=c(1,1))