Problem 3

library(KFAS)
library(Quandl)

This problem requires local level model for annual flow on river Nile. For this task we use in built Nile data of R package.

(a)

y <- datasets::Nile
str(y)
##  Time-Series [1:100] from 1871 to 1970: 1120 1160 963 1210 1160 1160 813 1230 1370 1140 ...
plot(y, xlab="", ylab=expression(""*y[t]*""), main="Measurements of the annual flow of the river Nile at Ashwan ")

Now let’s modify and re-estimate the model (without missing observations) so that it includes a dummy variable dam taking value 0 before the Ashwan Dam was built in 1898, and value 1 afterwards. And define the state-space local level model.

dam <- ifelse(index(Nile)>1898,1,0)
T <-length(y)
y.LLM <- SSModel(y ~
                     SSMtrend(1, Q=list(NA)) +
                     dam,
                     data=Nile, H=NA)
# maximum likelihood estimation of paramaters of Q and H (i.e. variance of the two inovations)
pars <- log( rep(var(y, na.rm=TRUE), 2) )
y.fit <- fitSSM( model=y.LLM, inits=pars, method='BFGS')
# run Kalman Filter and Smoother with estimated parameters
y.out <- KFS(y.fit$model)

# construct 90% confidence intervals for filtered state - shorter way, using predict function with missing n.ahead option
y.KF.CI <- predict(y.fit$model, interval="confidence", level=0.9, filtered=TRUE)
# construct 90% confidence intervals for smoothed state
y.KS.CI <- predict(y.fit$model, interval="confidence", level=0.9)

# replace the filtered state for first period by NA
y.KF.CI[1,] <- NA

(b)

# filtered
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
         col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","filtered","90% confidence interval"), col=c(1,4,4), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )

# smoothed
plot.ts( cbind(y, y.KS.CI, Nile), plot.type="single",
         col=c(1,2,2,2,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","smoothed","90% confidence interval"), col=c(1,2,2), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )