inst/BookEx/C14R2.R

In FRAPO: Financial Risk Modelling and Portfolio Optimisation with R

MCsim <- function(x0 = 0, n, theta0, theta1){
    ans <- vector()
    length(ans) <- n + 1
    ans[1] <- x0
    for(i in 2:(n + 1)){
        ans[i] <- theta0 + theta1 * ans[i-1] + rnorm(1)
        }
    ans
}
## Parameter settings
theta0 <- 2
theta1 <- 0.5
N <- 10000
x01 <- 14
x02 <- -10
## Markov Chains
mc1 <- MCsim(x0 = x01, n = N, theta0 = theta0, theta1 = theta1)
mc2 <- MCsim(x0 = x02, n = N, theta0 = theta0, theta1 = theta1)
## Progression of Markov Chains
plot(mc1[1:100], type = "l", ylim = range(cbind(mc1, mc2)),
     xlab = "", ylab = "X", main = "Progression of first-order Markov Chain")
lines(mc2[1:100], col = "red")
## Expected value of stationarity distribution and estimates
EfPop <- theta0 / (1 - theta1) 
m <- trunc(N / 2) ## burn-in
EfEst1 <- mean(mc1[-c(1:m)])
c(EfPop, EfEst1)
## Standard error of estimate for first MC
ar1Est <- ar(mc1, order.max = 1)$ar
se1 <- sqrt(var(mc1) / N * (1 + ar1Est) / (1 - ar1Est))
c(EfEst1 - 2 * se1, EfEst1, EfEst1 + 2 * se1)
## Variance of stationarity distribution and estimate
VarfPop <- 1 / (1 - theta1^2)
VarfEst1 <- 1 / (1 - ar1Est^2)
c(VarfPop, VarfEst1)