# Code examples for the simulation chapter of ADAfaEPoV
# See that chapter for context


############### Code Example 1 ##################
# Illustration of the rejection method

# Draw one sample by the rejection method from the target distribution
# Inputs: target probability density function (dtarget)
  # proposal p.d.f. (dproposal)
  # random number generator for proposal (rproposal)
  # scaling factor (M)
# Presumes: dtarget takes a number, and nothing else, and returns a number;
  # so does dproposal;
  # rproposal takes the number of variables desired, and nothing else
# Output: A single random number, drawn from the target distribution
rrejection.1 <- function(dtarget,dproposal,rproposal,M) {
  rejected <- TRUE
  while(rejected) {
    R <- rproposal(1)
    U <- runif(1)
    rejected <- (M*U*dproposal(R) < dtarget(R))
  }
  return(R)
}

# Draw multiple samples by the rejection method
  # Call rrejection.1() multiple times, using replicate()
# Inputs: number of samples to draw (n)
  # target probability density function (dtarget)
  # proposal p.d.f. (dproposal)
  # random number generator for proposal (rproposal)
  # scaling factor (M)
# Calls: rrejection.1
# Output: vector of numbers of length n
rrejection <- function(n,dtarget,dproposal,rproposal,M) {
  replicate(n,rrejection.1(dtarget,dproposal,rproposal,M))
}

# Comment on these functions:
# This is not "industrial strength" programming, because the arguments which
# are functions (dtarget, dproposal, rproposal) are only allowed to be
# functions of one argument.  If we wanted to draw from a Gaussian with mean 7
# and sd 1, we couldn't say dtarget=dnorm, or dtarget=dnorm(x,mean=7,sd=1).
# Instead, we'd have to first define a new one-argument function, as
### dnorm.7.1 <- function(x) {dnorm(x,mean=7,sd=1)}
# and then say dtarget=dnorm.  This is way too clunky.  R does have mechanisms
# for passing arguments along from one function to the other functions it 
# invokes, but they'd distract from getting across the idea of the method.



################# Figure 1 #################
# Illustration of the rejection method with a beta distribution
  # See the notes
M <- 3.3
curve(dbeta(x,5,10),from=0,to=1,ylim=c(0,M))
r <- runif(300,min=0,max=1)
u <- runif(300,min=0,max=1)
below <- which(M*u*dunif(r,min=0,max=1) <= dbeta(r,5,10))
points(r[below],M*u[below],pch="+")
points(r[-below],M*u[-below],pch="-")



###### Code Example 2 #####
# The Arnold cat map, and using it as a pseudo-random number generator

# Arnold cat map function
  # See text of lecture for definition
# Input: a vector v, presumed to be two-dimensional
# Output: the result of applying the map to the vector
arnold.map <- function(v) {
  theta <- v[1]
  phi <- v[2]
  theta.new <- (theta+phi)%%1
  phi.new <- (theta+2*phi)%%1
  return(c(theta.new,phi.new))
}

# Pseudo-random number generator based on the Arnold cat map
  # Iterates applying the map, starting with a "seed" value, and returns the
  # sequence of values for the first coordinate ("theta" in terms of the notes)
# Inputs: number of random variates desired (n)
  # Two-dimensional initial condition vector for the map (seed)
# Output: sequence of n numbers, (approximately) uniformly distributed between
  # 0 and 1
rarnold <- function(n,seed) {
  z <- vector(length=n)
  for (i in 1:n) {
    seed <- arnold.map(seed)
    z[i] <- seed[1]
  }
  return(z)
}



#### Figure 2 ####
# Take from Wikipedia, s.v. "Arnold's cat map".




##### Figure 3 #####
# The Arnold cat map as a random number generator, illustrated
par(mfrow=c(2,1))
# Get 1000 draws from this generator
z <- rarnold(1000,c(0.11124,0.42111))
# Histogram: it's flat
hist(z,probability=TRUE)
# Scatterplot of successive values --- no predictable pattern, which is exactly
# what we want
plot(z[-1],z[-1000],xlab=expression(Z[t]),ylab=expression(Z[t+1]))


##### Code Example 3 #####
# Sampling from the fitted joint density of population growth and investment
# rates
  # This is an example of sampling from a kernel density estimate by first
  # picking training points, and then sampling from a copy of the kernel
  # distribution for each point
  # It's hard-coded to use the popinv2 object --- can you modify it so that we
  # could use it with anything produced by npudens?
# Inputs: number of draws to make (n)
# Presumes: popinv2 exists and has the attributes of "dens", "ndim", "eval",
  # and "bw"
# Output: n*2 matrix of draws, with columns named "pop.growth.rate" and
  # "invest.rate"
rpopinv <- function(n) {
  n.train <- length(popinv2$dens) # Number of training points
  ndim <- popinv2$ndim # Number of dimensions
  points <- sample(1:n.train,size=n,replace=TRUE) # Which training points
    # are the basis for our sample?
  z <- matrix(0,nrow=n,ncol=ndim)
  for (i in 1:ndim) {
    coordinates <- popinv2$eval[points,i] # Copy the coordinates of our points
    # Sample from the kernel distribution centered at each point
      # Presumes a Gaussian kernel (correct here)
      # Uses the fact that the mean and sd arguments of rnorm() can be vectors
    z[,i] <- rnorm(n,coordinates,popinv2$bw[i])
  }
  colnames(z) <- c("pop.growth.rate","invest.rate")
  return(z)
}

#### Figure 4 ####
# A challenge problem from lecture #6, code omitted

#### Figure 5 #####
# Plotting the distribution of rate ratios
YoverX <- z[,"invest.rate"]/z[,"pop.growth.rate"]
plot(density(YoverX),xlab="Y/X",ylab="Probability density",main="")
rug(YoverX,side=1)


##### Figure 6 ####
# The previous figure revisited, using overplotting to show sampling variability
# The function one.plot() basically wraps together what we did for the previous
# figure into a single command, so that we can repeat it multiple times
plot(0,xlab="Y/X",ylab="Probability density",type="n",xlim=c(-1,10),
     ylim=c(0,0.3))
one.plot <- function() {
  zprime <- rpopinv(2000)
  YoverXprime <- zprime[,"invest.rate"]/zprime[,"pop.growth.rate"]
  density.prime <- density(YoverXprime)
  lines(density.prime,col="grey")
}
invisible(replicate(50,one.plot()))
  # Surrounding the replicate() in invisible() keeps it from automatically
  # echoing to the screen --- see help(invisible)



##### Figure 7 ####
# Old Faithful data plus OLS regression line
plot(geyser$duration,geyser$waiting,xlab="duration",ylab="waiting")
abline(fit.ols)


##### Code Example 4 #####
# Simulating from the OLS model fit to the geyser data
  # Holds the values of the independent variable duration fixed, then generates
  # new values of the dependent variable waiting from the linear model plus
  # Gaussian noise
# Inputs: None
# Presumes: fit.ols exists and has the right attributes
# Outputs: A data frame with 299 rows and 2 columns, named duration and
  # waiting
rgeyser <- function() {
  n <- nrow(geyser)
  sigma <- summary(fit.ols)$sigma
  new.waiting <- rnorm(n,mean=fitted(fit.ols),sd=sigma)
  new.geyser <- data.frame(duration=geyser$duration,waiting=new.waiting)
  return(new.geyser)
}


#### Figure 8 ####
# Comparing the empirical distribution of waiting times to that implied by the
# OLS geyser model
plot(density(geyser$waiting),xlab="waiting",main="",sub="")
lines(density(rgeyser()$waiting),lty=2)


#### Figure 9 ####
# Comparing the scatterplot of real data to the simulation of the OLS model
plot(geyser$duration,geyser$waiting,xlab="duration",ylab="waiting")
abline(fit.ols)
points(rgeyser(),pch=20,cex=0.5)


#### Figure 10 ####
# Conditional variance function, as estimated from OLS residuals, plus what the
# OLS model says the variance is
library(np)
var1 <- npreg(residuals(fit.ols)^2 ~ geyser$duration)
plot(var1,xlab="duration",ylab="Conditional variance of waiting time")
abline(h=summary(fit.ols)$sigma^2,lty=2)

##### Figure 11 ####
# The previous plot, plus applying the variance-function estimation method to
# realizations of the OLS method, to see whether the curve we found on the data
# is compatible with homoskedasticity
  # The plotting technique here is basically the same as in previous notes
  # using functions from np --- set up a grid of points at which to evaluate
  # the function, do so, plot the fitted values.  The function
  # one.var.func() wraps up drawing a realization from rgeyser, calculating
  # the variance function, and adding it to the plot.
plot(var1)
abline(h=summary(fit.ols)$sigma^2,lty=2)
duration.grid <- seq(from=min(geyser$duration),to=max(geyser$duration),
                     length.out=300)
one.var.func <- function() {
  fit <- lm(waiting ~ duration, data=rgeyser())
  var.func <- npreg(residuals(fit)^2 ~ geyser$duration)
  lines(duration.grid,predict(var.func,exdat=duration.grid),col="grey")
}
invisible(replicate(10,one.var.func())) # See comments to Figure 6



##### Reading in the S&P 500 data ####
sp <- read.csv("SPhistory.short.csv")
# We only want closing prices
sp <- sp[,7]
# The data are in reverse chronological order, which is weird for us
sp <- rev(sp)
# And in fact we only want log returns, i.e., difference in logged prices
sp <- diff(log(sp))


##### Code Example 5 #####
# Functions for applying method of moments estimation to the MA(1) model

# 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 (c), variance (v)
# Calls: ma.mm.objective
# Outputs: The return value of optim(), which begins with the parameters
  # minimizing the objective function (see help(optim) for more)
ma.mm.est <- function(c,v) {
  theta.0 <- c/v
  sigma2.0 <- v
  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)
  # real covariance (c)
  # real variance (v)
# Output: The Euclidean distance
ma.mm.objective <- function(params,c,v) {
  theta <- params[1]
  sigma2 <- params[2]
  c.pred <- theta*sigma2
  v.pred <- sigma2*(1+theta^2)
  return((c-c.pred)^2 + (v-v.pred)^2)
}

##### Code Example 6 #####
# 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),] # z[1,] should really be "z[0,]"
  return(x)
}


##### Code Example 7 #####
# Simulate an MA(1) model and find variance and covariance

# 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))
}


##### Figure 12 #####
# Plot covariance, variance, and their ratio as functions of theta in the
# MA(1) model
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)


#### Code Example 8 ####
# Functions for method of simulated moments estimation with an MA(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)
# Calls: ma.msm.objective()
# Output: the return value of optim(), including the estimated parameters,
  # see help(optim)
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)
# Calls: sim.cov(), sim.var()
# Output: the distance between the measured moments and the simulation
  # estimate of them
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)
}