/ #ARIMA #ARMA 

Time Series Cheatbook

1 Library

knitr::opts_chunk$set(echo = TRUE)

library("tidyverse") #ggplot and dplyr 
library("gridExtra") # combine plots
library("knitr") # for pdf
library("fpp2") #timeseries with autoplot and stuff
library("reshape2") #reshape the data
library("MASS") #StepAIC
library("astsa") #dataset oil and gas is present here
library("zoo") #dataset oil and gas is present here
library("forecast") # for forecasting time series
library("kernlab") #gausspr function
library("TSA") #Q3
library("tseries")

# The palette with black:
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
                "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

set.seed(12345)
options(scipen=999)

1.1 Generate two time series and use smoothing filter

set.seed(12345)

n = 100
x <- vector(length = n)
x2 <- vector(length = n)

x[1] <- 0
x[2] <- 0

#first series generation
for(i in 3:n){
  x[i] <- -0.8 * x[i-2] + rnorm(1,0,1)
}

#second series generation
for(i in 1:n){
  x2[i] <- cos(0.4*pi*i)
}

# smoothing filter function
smoothing_filter <- function(x){
v <- vector(length = length(x))
for(i in 5:length(x)){
  v[i] = 0.2 * (x[i] + x[i-1] + x[i-2] + x[i-3] + x[i-4])
}
return(v)
}

#generate smoothed series
smooth_x <- smoothing_filter(x)
smooth_x2 <- smoothing_filter(x2)

#adding everything to a dataframe
df <- cbind(x,x2,smooth_x,smooth_x2) %>% as.data.frame() %>% mutate(index=1:100)

ggplot(df, aes(x=index)) + 
  geom_line(aes(y=x, color="Original Time Series")) + 
  geom_line(aes(y=smooth_x, color="Smoothed Time Series")) + 
  ggtitle("Plot of 1st time series and its smoothed version") +
    scale_colour_manual("", breaks = c("Original Time Series", "Smoothed Time Series"),
                        values = c("#CC79A7", "#000000"))

ggplot(df, aes(x=index)) + 
  geom_line(aes(y=x2, color="Original Time Series")) + 
  geom_line(aes(y=smooth_x2, color="Smoothed Time Series")) + 
  ggtitle("Plot of 2ND time series and its smoothed version") +
    scale_colour_manual("", breaks = c("Original Time Series", "Smoothed Time Series"),
                        values = c("#CC79A7", "#000000"))

1.2 Casuality and Invertiblity.

Causality: ARMA(p,q) is causal iff roots \(\phi(z')=0\) are outside unit circle. eg: \(x_t = 0.4x_{t-1} + 0.3x_{t-2} + w_t\), roots are -> \(1-0.4B+0.3B^2\)

equation is: \(\phi(Z) = 1-4B+2B^2+0 B^3 +0 B^4 + B^5\)

z = c(1,-4,2,0,0,1)
polyroot(z)
## [1]  0.2936658+0.000000i -1.6793817+0.000000i  1.0000000-0.000000i
## [4]  0.1928579-1.410842i  0.1928579+1.410842i
any(Mod(polyroot(z))<=1)
## [1] TRUE

Invertible: ARMA(p,q) is causal iff roots \(\theta(z')=0\) are outside unit circle.

equation is: \(\theta(Z) = 1+3B^2+B^4-4B^6\)

z = c(1,0,3,0,1,0,-4)
polyroot(z)
## [1]  0.1375513+0.6735351i -0.1375513+0.6735351i -0.1375513-0.6735351i
## [4]  0.1375513-0.6735351i  1.0580446+0.0000000i -1.0580446+0.0000000i
any(Mod(polyroot(z))<=1)
## [1] TRUE

1.3 ACF and Theortical ACF

set.seed(54321)

series <- arima.sim(n = 100, list(ar = c(-3/4), ma = c(0,-1/9)))

acf(series)

acf(ARMAacf(ar = c(-3/4), ma =  c(0,-1/9), lag.max = 20))

2 Visualization, detrending and residual analysis of Rhine data.

The data set Rhine.csv contains monthly concentrations of total nitrogen in the Rhine River in the period 1989-2002.

2.1 ACF Plot

set.seed(12345)
rhine_data <- read.csv2("../../static/data/Rhine.csv")
rhine_data_ts <- ts(data = rhine_data$TotN_conc, 
                    start = c(1989,1), 
                    frequency = 12)

plot.ts(rhine_data_ts, main="Time Series of Nitrogen Concentration in Rhine")

lag.plot(rhine_data_ts,lags = 12)  

acf(rhine_data_ts)

#alternative
autoplot(rhine_data_ts) + ylab("Total Concentration") +xlab("Year") +
  ggtitle("Concentration of Nitrogen in Rhine vs. Year")

gglagplot(rhine_data_ts, lags = 1, set.lags = 1:12, color=FALSE) 

ggAcf(rhine_data_ts) + ggtitle("ACF for Total Nitrogen Concentration")

2.2 Detrending using linear regression

set.seed(12345)

rhine_lm_model <- lm(TotN_conc~Time, data=rhine_data)
plot(rhine_lm_model$residuals, type = 'l', main="Plot of Residual from the 
     linear model of Nitrogen Concentration")

acf(rhine_lm_model$residuals)

2.3 Detrending using kernel smoother

set.seed(12345)

model_smooth_lag_5 <- ksmooth(x = rhine_data$Time, y = rhine_data$TotN_conc, 
                              bandwidth=5)
model_smooth_lag_10 <- ksmooth(x = rhine_data$Time, y = rhine_data$TotN_conc, 
                               bandwidth=10)
model_smooth_lag_20 <- ksmooth(x = rhine_data$Time, y = rhine_data$TotN_conc, 
                               bandwidth=20)

model_smooth_lag_5_residual <- rhine_data$TotN_conc - model_smooth_lag_5$y
model_smooth_lag_10_residual <- rhine_data$TotN_conc - model_smooth_lag_10$y
model_smooth_lag_20_residual <- rhine_data$TotN_conc - model_smooth_lag_20$y

residual_df <- cbind(model_smooth_lag_5_residual, model_smooth_lag_10_residual, 
                     model_smooth_lag_20_residual, rhine_data$Time) %>% 
  as.data.frame()

colnames(residual_df) <- c("lag_5_residual", "lag_10_residual", 
                           "lag_20_residual", "Time")

ggplot(residual_df, aes(x=Time)) + 
  geom_line(aes(y=lag_5_residual, color="Lag 5 residual")) + 
  geom_line(aes(y=lag_10_residual, color="Lag 10 residual")) + 
  geom_line(aes(y=lag_20_residual, color="Lag 20 residual")) + 
  ggtitle("Residual vs. Time by Lag") +
    scale_colour_manual("", breaks = c("Lag 5 residual", "Lag 10 residual", 
                                       "Lag 20 residual"),
                        values = c("#CC79A7", "#000000", "#D55E00"))

acf(model_smooth_lag_5_residual)

acf(model_smooth_lag_10_residual)

acf(model_smooth_lag_20_residual)

2.4 Detrending using seasonal means model

set.seed(12345)

rhine_data_wide <- rhine_data
rhine_data_wide$dummy <- "1"
rhine_data_wide$Month <- paste0("Month_",rhine_data_wide$Month)
rhine_data_wide <- dcast(rhine_data_wide, 
                         formula = TotN_conc+Year+Time~Month, 
                         value.var = "dummy", fill = "0")

lm_model_month_lag <- lm(data=rhine_data_wide, 
                    TotN_conc~Time+Month_1+Month_2+Month_3+Month_4+Month_5+
                      Month_6+
                      Month_7+
                      Month_8+Month_9+Month_10+Month_11+Month_12)


plot(lm_model_month_lag$residuals, type = 'l', main="Plot of the Residuals vs. Time")

acf(lm_model_month_lag$residuals)

2.5 Model tuning using SetpAIC

set.seed(12345)

lm_model_month_lag_step <- stepAIC(lm_model_month_lag, 
                                   scope = list(upper = ~Time+Month_1+Month_2+
                                                        Month_3+Month_4+Month_5+
                                                        Month_6+Month_7+
                                                        Month_8+Month_9+Month_10+
                                                        Month_11+Month_12, 
                                                                  lower = ~1), 
                                   trace = TRUE, 
                                   direction="backward")
## Start:  AIC=-202.02
## TotN_conc ~ Time + Month_1 + Month_2 + Month_3 + Month_4 + Month_5 + 
##     Month_6 + Month_7 + Month_8 + Month_9 + Month_10 + Month_11 + 
##     Month_12
## 
## 
## Step:  AIC=-202.02
## TotN_conc ~ Time + Month_1 + Month_2 + Month_3 + Month_4 + Month_5 + 
##     Month_6 + Month_7 + Month_8 + Month_9 + Month_10 + Month_11
## 
##            Df Sum of Sq     RSS      AIC
## - Month_4   1     0.200  43.436 -203.249
## - Month_1   1     0.220  43.456 -203.170
## - Month_3   1     0.331  43.567 -202.743
## <none>                   43.237 -202.023
## - Month_2   1     1.440  44.677 -198.517
## - Month_11  1     2.305  45.541 -195.297
## - Month_5   1     3.274  46.511 -191.760
## - Month_10  1     3.401  46.637 -191.303
## - Month_9   1     7.853  51.089 -175.986
## - Month_6   1     8.215  51.452 -174.797
## - Month_7   1    14.321  57.557 -155.959
## - Month_8   1    16.488  59.725 -149.749
## - Time      1   118.387 161.624   17.499
## 
## Step:  AIC=-203.25
## TotN_conc ~ Time + Month_1 + Month_2 + Month_3 + Month_5 + Month_6 + 
##     Month_7 + Month_8 + Month_9 + Month_10 + Month_11
## 
##            Df Sum of Sq     RSS      AIC
## <none>                   43.436 -203.249
## - Month_1   1     0.640  44.077 -202.790
## - Month_3   1     0.851  44.288 -201.988
## - Month_11  1     2.235  45.671 -196.819
## - Month_2   1     2.706  46.142 -195.096
## - Month_5   1     3.355  46.791 -192.748
## - Month_10  1     3.502  46.938 -192.223
## - Month_9   1     8.868  52.304 -174.036
## - Month_6   1     9.317  52.753 -172.602
## - Month_7   1    16.912  60.348 -150.004
## - Month_8   1    19.636  63.072 -142.586
## - Time      1   118.194 161.630   15.506
colnames(lm_model_month_lag_step$model)
##  [1] "TotN_conc" "Time"      "Month_1"   "Month_2"   "Month_3"  
##  [6] "Month_5"   "Month_6"   "Month_7"   "Month_8"   "Month_9"  
## [11] "Month_10"  "Month_11"

3 Analysis of oil and gas time series

Weekly time series oil and gas present in the package astsa show the oil prices in dollars per barrel and gas prices in cents per dollar.

3.1 Checking Stationary

set.seed(12345)

data_oil <- astsa::oil
data_gas <- astsa::gas

ts.plot(data_oil, data_gas, gpars = list(col = c("black", "red")))

#alternative

autoplot(ts(cbind(data_oil, data_gas), start = 2000, frequency = 52)) + 
           ylab("Price of Oil and Gas") +xlab("Year") + 
           ggtitle("Price of Oil and Gas vs. Years")

3.2 Log transformation to fix stationary

set.seed(12345)

autoplot(ts(cbind(log(data_oil), log(data_gas)), start = 2000, frequency = 52)) + 
           ylab("Log of price of Oil and Gas") +xlab("Year") + 
           ggtitle("Log of price of Oil and Gas vs. Years")

3.3 Detrending using difference method

set.seed(12345)

autoplot(ts(diff(log(data_oil), differences = 1), start = 2000, frequency = 52)) + 
           ylab("Price of Oil") +xlab("Year") + 
           ggtitle("Price of Oil with Diff 1 vs. Years")

autoplot(ts(diff(log(data_gas), differences = 1), start = 2000, frequency = 52)) + 
           ylab("Price of Gas") +xlab("Year") + 
           ggtitle("Price of Gas with Diff 1 vs. Years")

ggAcf(diff(log(data_oil), differences = 1), data_oil)

ggAcf(diff(log(data_gas), differences = 1), data_gas)

3.4 Detrending using smoother

set.seed(12345)

oil_price_one_diff <- diff(log(data_oil), differences = 1) 
gas_price_one_diff <- diff(log(data_gas), differences = 1) 

df <- data.frame(oil_price_one_diff=as.matrix(oil_price_one_diff), 
           gas_price_one_diff = as.matrix(gas_price_one_diff), 
                      time=time(oil_price_one_diff))

df <- na.omit(df)

df$gas_price_one_diff = lag(df$gas_price_one_diff,1)
df$gas_price_two_diff = lag(df$gas_price_one_diff,2)
df$gas_price_three_diff = lag(df$gas_price_one_diff,3)


df <- na.omit(df)

df$smooth_one_week_lag <- ksmooth(x = df$oil_price_one_diff, 
                                  y = df$gas_price_one_diff, 
                                  bandwidth = 0.05, kernel = "normal")$y
df$smooth_two_week_lag <- ksmooth(x = df$oil_price_one_diff, 
                                  y = df$gas_price_two_diff, 
                                  bandwidth = 0.05, kernel = "normal")$y
df$smooth_three_week_lag <- ksmooth(x = df$oil_price_one_diff, 
                                    y = df$gas_price_three_diff, 
                                    bandwidth = 0.05, kernel = "normal")$y

df <- na.omit(df)

ggplot(data=df, aes(x=oil_price_one_diff, y = gas_price_one_diff)) + geom_point() +
  geom_line(aes(y= smooth_one_week_lag, color= "smooth_one_week_lag")) +
      scale_colour_manual("", breaks = c("smooth_one_week_lag"),
                        values = c("#CC79A7")) +
  ggtitle("Smoothed Plot of one week lag")

ggplot(data=df, aes(x=oil_price_one_diff, y = gas_price_two_diff)) + geom_point() +
    geom_line(aes(y= smooth_two_week_lag, color= "smooth_two_week_lag")) +
      scale_colour_manual("", breaks = c("smooth_two_week_lag"),
                        values = c("#56B4E9")) +
  ggtitle("Smoothed Plot of two week lag")

ggplot(data=df, aes(x=oil_price_one_diff, y = gas_price_three_diff)) + geom_point() +
    geom_line(aes(y= smooth_three_week_lag, color= "smooth_three_week_lag")) +
      scale_colour_manual("", breaks = c("smooth_three_week_lag"),
                        values = c("#D55E00")) +
  ggtitle("Smoothed Plot of three week lag")

3.5 Detrending using linear regression

set.seed(12345)

df$oil_price_two_diff = lag(df$oil_price_one_diff,2)
df$x_t_more_zero <- ifelse(df$oil_price_one_diff>0,"1","0")
lm_model_lag <- lm(data=df, formula = gas_price_one_diff~x_t_more_zero+
                     oil_price_one_diff+oil_price_two_diff)
summary(lm_model_lag)
## 
## Call:
## lm(formula = gas_price_one_diff ~ x_t_more_zero + oil_price_one_diff + 
##     oil_price_two_diff, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22855 -0.03191  0.00366  0.03692  0.37535 
## 
## Coefficients:
##                     Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        -0.002139   0.004588  -0.466     0.641    
## x_t_more_zero1      0.006265   0.007304   0.858     0.391    
## oil_price_one_diff  0.084851   0.077622   1.093     0.275    
## oil_price_two_diff  0.222373   0.050788   4.378 0.0000144 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05511 on 534 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.04571,    Adjusted R-squared:  0.04035 
## F-statistic: 8.526 on 3 and 534 DF,  p-value: 0.00001538
plot(lm_model_lag$residuals, type = 'l', main="Residual vs. Time")

ggAcf(lm_model_lag$residuals)

4 Linear Regressions on Necessarily Lagged Variables and Appropriate Correlation

4.1 Generate AR and using PACF

\(\phi_{33} = corr(X_{t-3}-f_p,X_t-f_p)\) where \(f_p=\sum_{j=1}^p \phi_j X_{\tau-j}\)

set.seed(12345)
x_t <- arima.sim(model = list(ar = c(0.8,-0.2,0.1)), n=1000)
actual_pacf_value <- pacf(x_t, plot = FALSE)$acf[3]
df <- data.frame(x_t = as.vector(x_t))
df$x_t_lag_1 <- lag(df$x_t,1)
df$x_t_lag_2 <- lag(df$x_t,2)
df$x_t_lag_3 <- lag(df$x_t,3)
df <- na.omit(df)

# building models and getting their residuals
model_1_res <- lm(x_t ~ x_t_lag_1 + x_t_lag_2, data = df)$residuals
model_2_res <- lm(x_t_lag_3 ~ x_t_lag_1 + x_t_lag_2, data = df)$residuals

# theortical pacf values
theotical_pacf_value <- cor(x = model_1_res, y = model_2_res, 
                            use = "na.or.complete")

cat("The theoretical and actual value of PACF are: ", theotical_pacf_value, 
    actual_pacf_value)
## The theoretical and actual value of PACF are:  0.1146076 0.1170643

4.2 Methods of Moments, Conditional Least Squares and Maximum Likelihood

set.seed(12345)
x_t <- arima.sim(model = list(ar = c(0.8,0.1)), n=100)

method_yule_walker <- ar(x_t, order = 2, method = "yule-walker", aic = FALSE)$ar
method_cls <- ar(x_t, order = 2, method = "ols", aic = FALSE)$ar
method_mle <- ar(x_t, order = 2, method = "mle", aic = FALSE)$ar

df <- data.frame(rbind(method_yule_walker, method_cls,method_mle))

kable(df, caption = "Comparison of parameters using different methods")
Table 4.1: Comparison of parameters using different methods
ar1 ar2
method_yule_walker 0.8029146 0.1037053
method_cls 0.8066782 0.1205352
method_mle 0.7968774 0.1189369
# Since varience is not given by ar we use arima function
ML_Model_CI = arima(x_t, order = c(2,0,0), method = "ML")
sigma = ML_Model_CI$var.coef[2, 2]
phi_2 = ML_Model_CI$coef[2]
CI = c(phi_2 - 1.96 * sigma, phi_2 + 1.96 * sigma)
CI
##        ar2        ar2 
## 0.09924714 0.13846032

5 ARIMA

5.1 Sample and Theoretical ACF and PACF

Now \(ARIMA(1,1,1)(1,1,1)_4\) can be written as \((1-\phi_1B)(1-B)(1-B^4)(1-\Phi_1 B^4)x_t = w_t (1+\theta B)(1+\Theta B^4)\) Similarly \(ARIMA(0,0,1)(0,0,1)_{12}\) can be written as \(x_t=w_t(1+\Theta B^{12})(1+\theta B)\) which can be simplified as \(x_t = w_t(1+\Theta B^{12}+ \theta B + \Theta \theta B^{13})\) given that \(\theta=0.3\) and \(\Theta =0.6\) we get \(x_t =w_t(1+0.3B+0.6B^{12}+0.18B^{13})\)

set.seed(12345)
x_t <- arima.sim(model = list(ma = c(0.3,rep(0,10),0.6,0.18)), n=200)

df <- data.frame(sample_acf = acf(x_t, plot = FALSE, lag.max = 14)$acf,
                 sample_pacf = pacf(x_t, plot = FALSE, lag.max = 14)$acf,
                 theortical_acf = ARMAacf(ma = c(0.3,rep(0,10),0.6,0.18), 
                                          pacf = FALSE, lag.max = 13),
                 theortical_pacf = ARMAacf(ma = c(0.3,rep(0,10),0.6,0.18), 
                                           pacf = TRUE, lag.max = 14))

df$index <- rownames(df)

plot1 <- ggplot(data=df, aes(x=index)) + 
  geom_col(aes(y=sample_acf)) + 
  ggtitle("Sample ACF")

plot2 <- ggplot(data=df, aes(x=index)) + 
  geom_col(aes(y=theortical_acf)) + 
  ggtitle("Theoretical ACF")

grid.arrange(plot1, plot2, ncol = 1)

plot3 <- ggplot(data=df, aes(x=index)) + 
  geom_col(aes(y=sample_pacf)) + 
  ggtitle("Sample PACF")

plot4 <- ggplot(data=df, aes(x=index)) + 
  geom_col(aes(y=theortical_pacf)) + 
  ggtitle("Theoretical PACF")

grid.arrange(plot3, plot4, ncol = 1)

5.2 Forecast and Prediction

set.seed(12345)
x_t <- arima.sim(model = list(ma = c(0.3,rep(0,10),0.6,0.18)), n=200)
fit_x_t <- arima(x_t, order = c(0,0,1), seasonal = list(order = c(0,0,1),
                                                        period = 12))
predicted_x_t <- predict(fit_x_t, n.ahead=30) 
predicted_x_t_upper_band <- predicted_x_t$pred + 1.96 * predicted_x_t$se
predicted_x_t_lower_band <- predicted_x_t$pred - 1.96 * predicted_x_t$se

#kernlab

df <- data.frame(y = x_t)
df$x <- as.numeric(rownames(df))
gausspr_model <- gausspr(x=df$x, y=df$y)
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
predicted_x_t_kernlab <- predict(gausspr_model, newdata=data.frame(x=201:230))

df3 <- data.frame(y = predicted_x_t_kernlab, x=201:230) 


df2 <- data.frame(predicted_x_t = predicted_x_t$pred, 
                  predicted_x_t_upper = predicted_x_t_upper_band, 
                  predicted_x_t_lower = predicted_x_t_lower_band,
                  x = 201:230)

ggplot() + 
  geom_line(data=df, aes(x=x, y=y, color="Actual y")) + 
    geom_line(data=df2, aes(x=x, y=predicted_x_t, color="Predicted y")) + 
      geom_line(data=df2, aes(x=x, y=predicted_x_t_upper, color="Upper band")) + 
        geom_line(data=df2, aes(x=x, y=predicted_x_t_lower, color="Lower band")) + 
      scale_colour_manual("", breaks = c("Actual y", "Predicted y", "Upper band", "Lower band"),
                        values = c("#000000", "#009E73", "#56B4E9", "#E69F00")) +
  ggtitle("Original vs. Predicted y with confidence bands")

ggplot() + 
  geom_line(data=df, aes(x=x, y=y, color="Actual y")) + 
    geom_line(data=df3, aes(x=x, y=y, color="Predicted y")) + 
        scale_colour_manual("", breaks = c("Actual y", "Predicted y"),
                        values = c("#000000","#56B4E9")) +
  ggtitle("Original vs. Predicted y using gausspr")

5.3 Prediction Band

x_t <- arima.sim(model = list(ar = c(0.7), ma=c(0.5)), n=50)
fit_x_t <- arima(x_t[1:40], order = c(1,0,1), include.mean = 0)

predicted_x_t <- predict(fit_x_t, n.ahead=10)
predicted_x_t_upper_band <- predicted_x_t$pred + 1.96 * predicted_x_t$se
predicted_x_t_lower_band <- predicted_x_t$pred - 1.96 * predicted_x_t$se

df <- data.frame(y = x_t[1:40], x=1:40) 
df2 <- data.frame(y = predicted_x_t$pred, 
                  upper_band=predicted_x_t_upper_band, 
                  lower_band=predicted_x_t_lower_band,
                  x = 41:50)

ggplot() + 
  geom_line(data=df, aes(x=x, y=y, color="Actual y")) + 
    geom_line(data=df2, aes(x=x, y=y, color="Predicted y")) + 
      geom_line(data=df2, aes(x=x, y=upper_band, color="Upper band")) + 
        geom_line(data=df2, aes(x=x, y=lower_band, color="Lower band")) + 
      scale_colour_manual("", breaks = c("Actual y", "Predicted y", "Upper band", "Lower band"),
                        values = c("#000000", "#009E73", "#56B4E9", "#E69F00")) +
  ggtitle("Original vs. Predicted y with confidence bands")  

5.4 Finding a Suitable ARIMA Model and EACF

set.seed(12345)

# visualization
autoplot(ts(oil, start = 2000, frequency = 52)) + 
           ylab("Price of Oil") +xlab("Year") + 
           ggtitle("Price of Oil vs. Years")

ggAcf(oil) + ggtitle("ACF for Oil")

ggAcf(diff(oil)) + ggtitle("ACF for Oil with one diff")

ggPacf(oil) + ggtitle("PACF for Oil")

ggPacf(diff(oil)) + ggtitle("PACF for Oil with one diff")

# with log
autoplot(ts(log(oil), start = 2000, frequency = 52)) + 
           ylab("Price of Oil in Log") +xlab("Year") + 
           ggtitle("Price of Log Oil vs. Years")

autoplot(ts(diff(log(oil), lag=1), start = 1948, frequency = 12)) + 
           ylab("# Log Oil") +xlab("Year") + 
           ggtitle("Price of log oil with one lags vs. Years")

ggAcf(log(oil)) + ggtitle("ACF for log Oil")

ggAcf(diff(log(oil))) + ggtitle("ACF for log Oil with one diff")

ggPacf(log(oil)) + ggtitle("PACF for log Oil")

ggPacf(diff(log(oil))) + ggtitle("PACF for log Oil with one diff")

# EACF
eacf(diff(log(oil)))
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x o o o o x o o o  o  o  o 
## 1 x o x o o o o x o o o  o  o  o 
## 2 x x x o o o o x o o o  o  o  o 
## 3 x x x o o o o x o o o  o  o  o 
## 4 x o x o o o o x o o o  o  o  o 
## 5 x x x o x o o x o o o  o  o  o 
## 6 o x x o x x o x o o o  o  o  x 
## 7 o x x x x x x x o x o  o  o  o

Analysis: ARIMA(0,1,1) or ARIMA(1,1,1) or ARIMA(0,1,3) according to EACF

#Suggested Models
modelA <- sarima(log(oil), 0,1,1)
## initial  value -3.058495 
## iter   2 value -3.068906
## iter   3 value -3.069474
## iter   4 value -3.069476
## iter   4 value -3.069476
## iter   4 value -3.069476
## final  value -3.069476 
## converged
## initial  value -3.069450 
## iter   2 value -3.069450
## iter   2 value -3.069450
## iter   2 value -3.069450
## final  value -3.069450 
## converged

modelB <- sarima(log(oil), 1,1,1)
## initial  value -3.057594 
## iter   2 value -3.061420
## iter   3 value -3.067360
## iter   4 value -3.067479
## iter   5 value -3.071834
## iter   6 value -3.074359
## iter   7 value -3.074843
## iter   8 value -3.076656
## iter   9 value -3.080467
## iter  10 value -3.081546
## iter  11 value -3.081603
## iter  12 value -3.081615
## iter  13 value -3.081642
## iter  14 value -3.081643
## iter  14 value -3.081643
## iter  14 value -3.081643
## final  value -3.081643 
## converged
## initial  value -3.082345 
## iter   2 value -3.082345
## iter   3 value -3.082346
## iter   4 value -3.082346
## iter   5 value -3.082346
## iter   5 value -3.082346
## iter   5 value -3.082346
## final  value -3.082346 
## converged

modelC <- sarima(log(oil), 0,1,3)
## initial  value -3.058495 
## iter   2 value -3.086110
## iter   3 value -3.086980
## iter   4 value -3.087501
## iter   5 value -3.087521
## iter   6 value -3.087521
## iter   7 value -3.087522
## iter   8 value -3.087522
## iter   9 value -3.087522
## iter   9 value -3.087522
## iter   9 value -3.087522
## final  value -3.087522 
## converged
## initial  value -3.087448 
## iter   2 value -3.087448
## iter   3 value -3.087449
## iter   3 value -3.087449
## iter   3 value -3.087449
## final  value -3.087449 
## converged

#ADF test
adf.test(modelA$fit$residuals)
## Warning in adf.test(modelA$fit$residuals): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  modelA$fit$residuals
## Dickey-Fuller = -6.5298, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
adf.test(modelB$fit$residuals)
## Warning in adf.test(modelB$fit$residuals): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  modelB$fit$residuals
## Dickey-Fuller = -6.461, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
adf.test(modelC$fit$residuals)
## Warning in adf.test(modelC$fit$residuals): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  modelC$fit$residuals
## Dickey-Fuller = -6.7187, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
#Redundancy check
summary(modelA$fit)
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ma1  constant
##       0.1701    0.0018
## s.e.  0.0499    0.0023
## 
## sigma^2 estimated as 0.002157:  log likelihood = 897.88,  aic = -1789.76
## 
## Training set error measures:
## Warning in trainingaccuracy(f, test, d, D): test elements must be within
## sample
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN
summary(modelB$fit)
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1     ma1  constant
##       -0.5264  0.7146    0.0018
## s.e.   0.0871  0.0683    0.0022
## 
## sigma^2 estimated as 0.002102:  log likelihood = 904.89,  aic = -1801.79
## 
## Training set error measures:
## Warning in trainingaccuracy(f, test, d, D): test elements must be within
## sample
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN
summary(modelC$fit)
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ma1      ma2     ma3  constant
##       0.1688  -0.0900  0.1447    0.0017
## s.e.  0.0424   0.0425  0.0430    0.0024
## 
## sigma^2 estimated as 0.00208:  log likelihood = 907.67,  aic = -1805.34
## 
## Training set error measures:
## Warning in trainingaccuracy(f, test, d, D): test elements must be within
## sample
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN
#BIC
BIC(modelA$fit)
## [1] -1776.859
BIC(modelB$fit)
## [1] -1784.592
BIC(modelC$fit)
## [1] -1783.844
#AIC
AIC(modelA$fit)
## [1] -1789.756
AIC(modelB$fit)
## [1] -1801.787
AIC(modelC$fit)
## [1] -1805.339
#Model C is the best

According to AIC and BIC the ModelC (ARIMA 0,1,3) is the best

Model equation is \(\Delta x_t=w_t+0.1688w_{t-1}-0.0900w_{t-2}+0.1447w_{t-3}\)

#Forecasting
sarima.for(log(oil), 0,1,3, n.ahead = 20)

## $pred
## Time Series:
## Start = c(2010, 26) 
## End = c(2010, 45) 
## Frequency = 52 
##  [1] 4.222141 4.222731 4.212938 4.214647 4.216356 4.218066 4.219775
##  [8] 4.221485 4.223194 4.224904 4.226613 4.228323 4.230032 4.231741
## [15] 4.233451 4.235160 4.236870 4.238579 4.240289 4.241998
## 
## $se
## Time Series:
## Start = c(2010, 26) 
## End = c(2010, 45) 
## Frequency = 52 
##  [1] 0.04561249 0.07016150 0.08569792 0.10226755 0.11650396 0.12918085
##  [7] 0.14072033 0.15138273 0.16134203 0.17072132 0.17961149 0.18808192
## [13] 0.19618697 0.20397021 0.21146718 0.21870731 0.22571532 0.23251220
## [19] 0.23911597 0.24554218

6 State Space Models

In table 1 a script for generation of data from simulation of the following state space model and implementation of the Kalman filter on the data is given.

\[ Z_t = A_{t-1} Z_{t-1} + e_t\] \[x_t = C_t z_t + \nu_t\]

\[\nu_t \sim N(0,R_t)\] \[e_t \sim N(0,Q_t)\]

\[ z_k = z_{k-1} + N(0,Q_t)\] \[x_k = z_{k} + N(0, R_t)\]

# generate  dataset
set.seed(1)
num = 50
w = rnorm(num+1,0,1)
v = rnorm(num ,0,1)
mu = cumsum(w) # state: mu[0], mu[1],..., mu[50]
y = mu[-1] + v # obs: y[1],..., y[50]
# filter  and  smooth (Ksmooth0 does  both)
ks = Ksmooth0(num , y, A=1, mu0=0, Sigma0=1, Phi=1, cQ=1, cR=1)
# start  figurepar(mfrow=c(3,1))
Time = 1:num
plot(Time , mu[-1], main="Predict", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xp)
lines(ks$xp+2*sqrt(ks$Pp), lty=2, col=4)
lines(ks$xp -2*sqrt(ks$Pp), lty=2, col=4)

plot(Time , mu[-1], main="Filter", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xf)
lines(ks$xf+2*sqrt(ks$Pf), lty=2, col=4)
lines(ks$xf -2*sqrt(ks$Pf), lty=2, col=4)

plot(Time , mu[-1], main="Smooth", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xs)
lines(ks$xs+2*sqrt(ks$Ps), lty=2, col=4)
lines(ks$xs-2*sqrt(ks$Ps), lty=2, col=4)

mu[1]
## [1] -0.6264538
ks$x0n
##            [,1]
## [1,] -0.3241541
sqrt(ks$P0n) # initial  value  info
##           [,1]
## [1,] 0.7861514
# filering with moving average 5
plot(Time , mu[-1], ylim=c(-5,10), main="Moving average smoothing with order 5")
lines(Time ,y,col="green")
lines(ks$xf)
lines(ks$xf+2*sqrt(ks$Pf), lty=2, col=4)
lines(ks$xf -2*sqrt(ks$Pf), lty=2, col=4)
lines(ma(y, order=5, ), col=6)

Analysis: We find that the moving average smoothing function with order 5 is the worst fit since its losing all the variability that is captured by our kalman filter.

# filter  and  smooth (Ksmooth0 does  both)
ks = Ksmooth0(num , y, A=1, mu0=0, Sigma0=1, Phi=1, cQ=10, cR=0.1)
# start  figurepar(mfrow=c(3,1))
Time = 1:num
plot(Time , mu[-1], main="Predict", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xp)
lines(ks$xp+2*sqrt(ks$Pp), lty=2, col=4)
lines(ks$xp -2*sqrt(ks$Pp), lty=2, col=4)

plot(Time , mu[-1], main="Filter", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xf)
lines(ks$xf+2*sqrt(ks$Pf), lty=2, col=4)
lines(ks$xf -2*sqrt(ks$Pf), lty=2, col=4)

Analysis: We find that the filtering output resembles the true value much more than the previousily run values.

# filter  and  smooth (Ksmooth0 does  both)
ks = Ksmooth0(num , y, A=1, mu0=0, Sigma0=1, Phi=1, cQ=0.1, cR=10)
# start  figurepar(mfrow=c(3,1))
Time = 1:num
plot(Time , mu[-1], main="Predict", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xp)
lines(ks$xp+2*sqrt(ks$Pp), lty=2, col=4)
lines(ks$xp -2*sqrt(ks$Pp), lty=2, col=4)

plot(Time , mu[-1], main="Filter", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(ks$xf)
lines(ks$xf+2*sqrt(ks$Pf), lty=2, col=4)
lines(ks$xf -2*sqrt(ks$Pf), lty=2, col=4)

Analysis: We find that the filtering output is far off the true value and hardly moving in terms of the predicted values, this is because kalman gain assumes that the uncertainity in measurement is high and it largely sticks with the mean of the series.

6.1 Kalman Filter Code

#Times: Number of time steps
#A: How the mean of z_t will be affected by z_(t-1) Used in transition model
#C: Scale the mean (z_t) in the emission model
#Q: Covariate matrix in the emission model
#R: Covariate matrix in the transition model 
#y: All observations
#mu0: first mean
#sigma0: first varience
my_kalman <- function(Times, y, A, mu0, Sigma0, C, Q, R) {
  mu <- rep(1,Times)
  sigma <- rep(1,Times)
  my_unweighted <- rep(1,Times)
  sigma_unweighted <- rep(1,Times)
  kalman_gain <- rep(1,Times)
  
  # Our best guess is that my_1 is mu0 else it could be the first observation
  mu[1] <- mu0
  
  # We don't know what sigma_1 is, so we choose 1 else it would be provided
  sigma[1] <- Sigma0
  
  for (t in 2:Times) {
  # Calculate the unweighted prediction of the mean
  my_unweighted[t] <- A[t]%*%mu[t - 1]
    
  # Calculate the unweighted prediction of the covariate matrix
  sigma_unweighted[t] <- A[t]%*%sigma[t - 1]%*%t(A[t]) + Q[t]
    
  # Kalman gain Used to weight between our unweighted prediction and the obs
  kalman_gain[t] <- sigma_unweighted[t]%*%t(C[t]) %*% 
    solve(C[t]%*%sigma_unweighted[t]%*%t(C[t] + R[t]))
  
  # Calculate the weighted mean, thus our prediction of the hidden state
  mu[t] <- my_unweighted[t] + kalman_gain[t]%*%(y[t] - C[t]%*%my_unweighted[t])
    
    # Calculate the weighted covariance matrix, thus our prediction of the predition error
    sigma[t] <- (1 - kalman_gain[t]%*%C[t])%*%sigma_unweighted[t]
  }
  
  return (list(mu = mu, sigma = sigma))
}

6.1.1 Alternative version

### MODIFIED, DOES WORK

kalman_filter = function(num, data, A, C, Q, R, m0, P0) {
  
  if (length(A) == 1)
    A = as.list(rep(A, num+1))

  if (length(C) == 1)
    C = as.list(rep(C, num+1))
  
  if (length(Q) == 1)
    Q = as.list(rep(Q, num+1))
  
  if (length(R) == 1)
    R = as.list(rep(R, num+1))
  
  # Init
  m = list()
  P = list()
  K = list()
  
  # Setup
  # Note that the formula with the inverse has been changed, as otherwise the
  # dimensions have to be determined first.
  m[[1]] = m0
  P[[1]] = P0
  
  for (t in 2:(num)) {
    
    # Prediction Step
    m[[t-1]] = A[[t-1]] %*% m[[t-1]]
    P[[t-1]] = A[[t-1]] %*% P[[t-1]] %*% t(A[[t-1]]) + Q[[t]]
    
    # Observation Update Step
    K[[t]] = P[[t-1]] %*% t(C[[t]]) %*% solve(C[[t]] %*% P[[t-1]] 
                                              %*% t(C[[t]]) + R[[t]])
    m[[t]] = m[[t-1]] + K[[t]] %*% (data[[t]] - C[[t]] %*% m[[t-1]])
    P[[t]] = P[[t-1]] - K[[t]] %*% C[[t]] %*% P[[t-1]]
    #P[[t]] = (diag(ncol(K[[t]])) - K[[t]] %*% C[[t]]) %*% P[[t-1]]
    
  }
  
  return(list(m = unlist(m), P=unlist(P), K=unlist(K)))
}

# generate dataset
set.seed(1)
num = 50
w = rnorm(num+1, 0, 1)
v = rnorm(num, 0, 1)
mu = cumsum(w) # state: mu[0], mu[1],..., mu[50]
y = mu[-1] + v # obs: y[1],..., y[50]

# filter  and  smooth (Ksmooth0 does both)
ks = kalman_filter(num, y, A=1, m0=0, P0=1, C=1, Q=1, R=1)
# start figurepar(mfrow=c(3,1))
Time = 1:num

plot(Time, mu[-1], main="Predict (custom Kalman Filter)", ylim=c(-5, 10))
lines(Time, y, col="green")
lines(ks$m)
lines(ks$m+2*sqrt(ks$P), lty=2, col=4)
lines(ks$m-2*sqrt(ks$P), lty=2, col=4)

# generate  dataset
set.seed(1)
num = 50
w = rnorm(num+1,0,1)
v = rnorm(num ,0,1)
mu = cumsum(w) # state: mu[0], mu[1],..., mu[50]
y = mu[-1] + v # obs: y[1],..., y[50]



my_ks = my_kalman(Times=num , y=y, A=rep(1,num), mu0=0, Sigma0=1, 
                  C=rep(1,num), Q=rep(1,num), R=rep(1,num))



plot(Time , mu[-1], main="Predict", ylim=c(-5,10))
lines(Time ,y,col="green")
lines(my_ks$mu)
lines(my_ks$mu+2*sqrt(my_ks$sigma), lty=2, col=4)
lines(my_ks$mu -2*sqrt(my_ks$sigma), lty=2, col=4)

Analysis: Kalman gain is given by \(K=\frac{P_{k} H^T_{k}}{P_{k} H^T_{k} + R_k}\) where you will realize that the relative magnitudes of matrices \(R_k\) and \(P_k\) control a relation between the filter’s use of predicted state estimate \(z_t\) and measurement \(x_t\).

When \(R_k\) tends to zero then \(x_t = x_{t-1} + K (y_t - H_k)\) suggests that when the magnitude of R is small, meaning that the measurements are accurate, the state estimate depends mostly on the measurements.

When the state is known accurately, then numerator is small compared to R, and the filter mostly ignores the measurements relying instead on the prediction derived from the previous state.