## ----include=FALSE------------------------------------------------------------ ##### _Advanced Data Analysis from an Elementary Point of View_ ##### # R code for the chapter "Simulation-Based Inference" # Please do not re-distribute or use without attribution # http://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/ ## ----s-and-p-stock-example-for-simulation-inference--------------------------- sp <- pdfetch_YAHOO("SPY", fields="adjclose", from=as.Date("1993-02-09"), to=as.Date("2018-02-09")) # We only want log returns, i.e., difference in logged prices sp <- diff(log(sp)) # Drop the initial NA, to avoid difficulties later sp <- sp[-1] ## ----------------------------------------------------------------------------- # Method of moments estimator for MA(1) # Makes a crude initial guess about the parameters, then minimizes the # method-of-moments objective function, using the built-in optimization # function optim() # Inputs: covariance between successive observations (c) # variance of observations (v) # Output: optim() object giving best-fitting parameters, objective function, # diagnostics, etc. # Calls: ma.mm.objective() # Presumes: time-series is stationary, generated by MA(1) ma.mm.est <- function(c,v) { # Crude initial estimates of MA(1) parameters from the moments theta.0 <- c/v sigma2.0 <- v # Adjust the parameters to minimizing the method-of-moments objective fit <- optim(par=c(theta.0,sigma2.0), fn=ma.mm.objective, c=c, v=v) return(fit) } # Calculate the method of moments objective function for MA(1) # Euclidean distance between actual covariance and variance, and that # implied by the MA(1) model # Inputs: vector of MA(1) parameters (params) # covariance between successive observations (c) # variance of the observations (v) # Output: The squared Euclidean distance between real and implied moments # Presumes: params is a length-two vector, with theta first ma.mm.objective <- function(params,c,v) { # Extract the parameters from params --- presumes a fixed oder theta <- params[1] sigma2 <- params[2] # Calculate the covariance and variance implied by the parameters c.pred <- theta*sigma2 v.pred <- sigma2*(1+theta^2) # Return the sum-of-squared-differences return((c-c.pred)^2 + (v-v.pred)^2) } ## ----------------------------------------------------------------------------- # Simulate an MA(1) model # Generate z values, then add them up with appropriate weights to get the x # values, which are returned # Set up to allow for doing multiple parallel simulation runs # In principle, arima.sim() in the stats package can be used instead, but # that would be a mere black box # Inputs: length of the simulation run (n) # MA(1) persistence parameter (theta) # Variance of z terms (sigma2) - note not their standard deviation # Number of independent simulation runs to do (s) # Outputs: an n by s matrix where each column is an independent MA(1) series rma <- function(n,theta,sigma2,s=1) { z <- replicate(s, rnorm(n=n+1, mean=0, sd=sqrt(sigma2))) # n+1 because x[1,] needs two z values # replicate(s, foo): if each run of foo gives a vector, replicate makes # a matrix with s columns x <- z[-1,] + theta*z[-(n+1),] return(x) } ## ----------------------------------------------------------------------------- # Find the variance of an MA(1) by simulation # Inputs: length of simulation (n) # MA(1) memory parameter (theta) # Noise variance sigma2 # Number of independent simulation paths to average over (s) # Output: the variance sim.var <- function(n,theta,sigma2,s=1) { vars <- apply(rma(n,theta,sigma2,s),2,var) # Get variance of each column # running var on a matrix also calculates covariances across columns, # which we don't want return(mean(vars)) } # Find the covariance of an MA(1) by simulation # Inputs: length of simulation (n) # MA(1) memory parameter (theta) # Noise variance sigma2 # Number of independent simulation paths to average over (s) # Output: the covariance sim.cov <- function(n,theta,sigma2,s=1) { x <- rma(n,theta,sigma2,s) covs <- colMeans(x[-1,]*x[-n,]) return(mean(covs)) } ## ----cov-var-and-cor-vs-theta,echo=FALSE-------------------------------------- par(mfrow=c(2,2)) theta.grid <- seq(from=-1,to=1,length.out=300) cov.grid <- sapply(theta.grid,sim.cov,sigma2=1,n=length(sp),s=10) plot(theta.grid,cov.grid,xlab=expression(theta),ylab="Covariance") abline(0,1,col="grey",lwd=3) var.grid <- sapply(theta.grid,sim.var,sigma2=1,n=length(sp),s=10) plot(theta.grid,var.grid,xlab=expression(theta),ylab="Variance") curve((1+x^2),col="grey",lwd=3,add=TRUE) plot(theta.grid,cov.grid/var.grid,xlab=expression(theta), ylab="Ratio of covariance to variance") curve(x/(1+x^2),col="grey",lwd=3,add=TRUE) par(mfrow=c(1,1)) ## ----------------------------------------------------------------------------- # Method of simulated moments estimator for the MA(1) model # Start with a crude parameter guess, then minimize the objective function # using simulation to approximate moments, with minimization via # optim() # Inputs: sample covariance (c) # sample variance (v) # length of simulation paths (n) # number of simulation paths per parameter value (s) # Output: optim() object giving best-fitting parameters, objective function, # diagnostics, etc. # Calls: ma.msm.objective() ma.msm.est <- function(c,v,n,s) { theta.0 <- c/v sigma2.0 <- v fit <- optim(par=c(theta.0,sigma2.0),fn=ma.msm.objective,c=c,v=v,n=n,s=s) return(fit) } # Method of simulated moments objective function # Get implied values of the moments by simulation and see how far they are # from the data # Inputs: length-two vector of parameters (params) # sample covariance (c) # sample variance (v) # length of simulation paths (n) # number of simulation paths per parameter value (s) # Output: the squared distance between the measured moments and the simulation # estimate of them # Calls: sim.cov(), sim.var() ma.msm.objective <- function(params,c,v,n,s) { theta <- params[1] sigma2 <- params[2] c.pred <- sim.cov(n,theta,sigma2,s) v.pred <- sim.var(n,theta,sigma2,s) return((c-c.pred)^2 + (v-v.pred)^2) }