4.1 Filter the group to be aforecasted (S01..S06)

Lets create a variable just with values, and exclude last 140 rows which are the ones to be forecaste and are NA

varsize <- length(S01[,"Var01"]) - 140
VarA <- ts(S01[1:varsize,"Var01"])
VarB <- ts(S01[1:varsize,"Var02"])

4.2 Impute missing values

We will use interpolation to determine value missing.

VarA <- na.interp(VarA)
sum(is.na(VarA))

## [1] 0

VarB <- na.interp(VarB)
sum(is.na(VarB))

## [1] 0

4.3 Let’s adjust for outliers

VarA %>% autoplot()

VarB %>% autoplot()

There are 3 clear outliers in Var02, I think we could smooth the series too

outlier <- (which.max(VarB))
outlier

## [1] 580

We will replace the value with the MEAN of the values before and after

VarB[outlier] <- (VarB[outlier-1] + VarB[outlier+1])/2

We will do again for second outlier

outlier <- (which.max(VarB))
VarB[outlier] <- (VarB[outlier-1] + VarB[outlier+1])/2

We will do again for third outlier

outlier <- (which.max(VarB))
VarB[outlier] <- (VarB[outlier-1] + VarB[outlier+1])/2

Plot again

VarA %>% autoplot()

VarB %>% autoplot()

4.4 Let’s smooth the series

We defined our own “smooth function” as it later would give us to flexibility how we smooth things. In this case we use a simple Moving Average

MySmooth <- function(series, size) {

# Because MA creates NA's at both ends we will save the tails so we can re-introduce them afterwards
  lengthS <- length(series)
  tail <- (size-1)/2
  tail1 <- series[1:tail]
  tail2 <- series[lengthS-tail+1:lengthS]

# Now that tails are saved, les use MA to smooth the series
  tempseries <- ma(series,size, centre = TRUE)

    # Lets add the tails back

  t <- length(tail1)
  
  for (i in 1:t) {
    tempseries[i] <- tail1[i]
    tempseries[lengthS-i+1] <- tail2[t-i+1]
    }
  
  return(tempseries)
}

We will use MA3 for a little touch

VarA <- MySmooth(VarA,3)
VarB <- MySmooth(VarB,3)

4.5 Transform the data for non-normal distributions

VarA %>% density() %>% plot()

VarB %>% density() %>% plot()

Both of them have skewness

let’s transform them using BoxCox

lambdaA <- BoxCox.lambda(VarA)
lambdaB <- BoxCox.lambda(VarB)

VarA <- BoxCox(VarA, lambdaA)
VarB <- BoxCox(VarB, lambdaB)

plot again

VarA %>% density() %>% plot()

VarB %>% density() %>% plot()

6.1 FIt model using ARIMA

VarA

fitA <- auto.arima(VarA)
checkresiduals(fitA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,5) with drift
## Q* = 5.1956, df = 3, p-value = 0.158
## 
## Model df: 7.   Total lags used: 10

fitA

## Series: VarA 
## ARIMA(1,1,5) with drift 
## 
## Coefficients:
##          ar1     ma1     ma2      ma3      ma4     ma5  drift
##       0.9136  0.1579  0.0188  -0.8912  -0.0772  0.0576  1e-03
## s.e.  0.2236  0.2251  0.2410   0.2227   0.0254  0.0292  6e-04
## 
## sigma^2 = 6.312e-05:  log likelihood = 5538.14
## AIC=-11060.29   AICc=-11060.2   BIC=-11017.16

fitA$aic

## [1] -11060.29

Will accept the recommended model (1,1,5) with drift

fcastA <- fitA %>% forecast(140) 
fcastA %>% autoplot()

VarB

fitB <- auto.arima(VarB)
checkresiduals(fitB)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,4)(0,0,1)[11]
## Q* = 19.974, df = 14, p-value = 0.131
## 
## Model df: 8.   Total lags used: 22

fitB

## Series: VarB 
## ARIMA(3,1,4)(0,0,1)[11] 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2      ma3     ma4     sma1
##       1.2261  -0.3007  0.0091  -0.8842  0.0458  -0.7573  0.6051  -0.0510
## s.e.  0.0818   0.0609  0.0385   0.0778  0.0333   0.0263  0.0602   0.0272
## 
## sigma^2 = 0.4879:  log likelihood = -1715.98
## AIC=3449.97   AICc=3450.08   BIC=3498.49

fitB$aic

## [1] 3449.969

We will keep the selected model (3,1,4)(0,0,1)

#fitB <- Arima(VarB, order=c(2,1,3), seasonal = #c(1,0,0),include.drift = TRUE)
#checkresiduals(fitB)
#fitB$aic

6.2 Forecast values with chosen ARIMA model

fcastB <- fitB %>% forecast(140)
fcastB %>% autoplot()

6.3 Refit using Exogenous variable

We will refit Var01 with VarB as exogenous predictor

fitAx <- auto.arima(VarA, xreg = VarB)
checkresiduals(fitAx)

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(1,1,5) errors
## Q* = 4.6381, df = 3, p-value = 0.2003
## 
## Model df: 8.   Total lags used: 11

fitAx

## Series: VarA 
## Regression with ARIMA(1,1,5) errors 
## 
## Coefficients:
##          ar1     ma1     ma2      ma3      ma4     ma5  drift   xreg
##       0.9145  0.1664  0.0196  -0.8831  -0.0908  0.0523  1e-03  8e-04
## s.e.  0.3441  0.3456  0.3731   0.3441   0.0270  0.0341  6e-04  2e-04
## 
## sigma^2 = 6.243e-05:  log likelihood = 5547.88
## AIC=-11077.77   AICc=-11077.65   BIC=-11029.25

fitAx$aic

## [1] -11077.77

Better AIC with exogenous variable

Let’s get the forecasted numbers

fcastA <- forecast(fitAx, h= 140, xreg = fcastB$mean)
autoplot(fcastA)

fcastA %>% autoplot()

6.4 Transform back (Inverse BoxCox)

fcastVarA <- as.numeric(fcastA$mean)
fcastVarA <- InvBoxCox(fcastVarA, lambda = lambdaA)

fcastVarB <- as.numeric(fcastB$mean)
fcastVarB <- InvBoxCox(fcastVarB, lambda = lambdaB)