1. Regressing TS to TS

Sometimes, you want to model your TS lniarly dependaent on other time series.

Sales depends on Global Economy

Utility Demand depends on Temparature

Level of toxin depends on local industory production

In that case, your Model is that the time series of interest \(X_t\), has linear relstionship with your independent series \(B_t\), plus stationary noise, \[ X_t \hspace{3mm} = \hspace{3mm} a + b B_t + Y_t \hspace{10mm} Y_t \sim ARMA \] where \(x_t\) contains the “explanatory” TS.




Sales Data and Indicator

  #- Indicator and Sales -
  D  = read.table("https://nmimoto.github.io/datasets/ls2.csv", header=T)
  head(D)   # D has two columns
##   Index Sales
## 1 10.01 200.1
## 2 10.07 199.5
## 3 10.32 199.4
## 4  9.75 198.9
## 5 10.33 199.0
## 6 10.13 200.2
  Index = ts(D[,"Index"], start=c(1,1), freq=1)
  Sales = ts(D[,"Sales"], start=c(1,1), freq=1)

  #- plot of Sales and index
  plot(Index, type='o')

  plot(Sales, type='o')

  #- Regress B on A
  Reg1 = lm(Sales~Index)
  summary(Reg1)
## 
## Call:
## lm(formula = Sales ~ Index)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.7736  -4.2953   0.2078   5.1369  13.5720 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.8812     5.3308   5.793 4.01e-08 ***
## Index        16.8060     0.4476  37.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.644 on 148 degrees of freedom
## Multiple R-squared:  0.905,  Adjusted R-squared:  0.9043 
## F-statistic:  1410 on 1 and 148 DF,  p-value: < 2.2e-16
  plot(Sales, type='o', main='Sales (Blk) vs Scaled Index (Red)')
  lines(Reg1$coef[1] + Index*Reg1$coef[2], col="red")

  plot(Reg1$residuals, type="o", main="Difference between Blk and Red")
  abline(h=0)

  layout(matrix(1:2, 1, 2))
  acf(Reg1$residuals)
  pacf(Reg1$residuals)

  #- Model the residual

We must worry about what’s called spurious regression

This is also called Cointegration Problem

Residuals out of TS-TS regression must be stationary

2. Using COR for Two Stationary TS


Two Independent TS

What happens when two stationary time series that has nothing to do with each other are regressed?

Regression estimates comes directry from sample correlation. What does correlation say?

Simulate with two independent ARMA(1,1)

  #-- X1 and X2 are both ARMA(1,1), but independent  --
  X1  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24) ))
  X2  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24) ))

  #-- They have nothing to do with each other, and correlation says so --
  plot(X1, type="l", main=paste("cor=", round(cor(X1,X2), 4) ))
  lines(X2, col="red")

  cor(X1,X2)
## [1] 0.01241742

To see the overall behavior, repleat 1000 times and plot histogram of correlations

  #-- Repeat above in lopp to see overall behavior --
  n = 250
  COR1 = 0
  for (i in 1:1000) {

    X1  = arima.sim(n = n, list(ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = n, list(ar = c(0.8), ma = c(-.24) ))

    COR1[i] = cor(X1,X2)

  }

  hist(COR1, main="Cor of 1000 pairs of Ind St. TS", xlim=c(-1,1) )

  quantile(COR1, c(.05,.95))   #- Simulated empirical CI for Correlation of 0.
##         5%        95% 
## -0.1950126  0.1968019


Dependent TS

Repeat above, but make two ARMA dependent by using correlated \(\epsilon_t\)

  #- generate 2-D e from bivariate normal
  library(MASS)
  e = mvrnorm(n = 250, mu=c(0,0), Sigma=matrix(c(1,.7,.7,1), 2,2))
  head(e)
##             [,1]       [,2]
## [1,] -0.34752539 -0.3164672
## [2,] -0.24658567 -0.8746511
## [3,]  1.39731600  1.5123602
## [4,]  0.00651515  0.3116654
## [5,]  1.68671607  1.5654501
## [6,] -0.15292334 -0.1188732
  plot(e[,1], e[,2])

  #-- Generate X1 and X2 with correlated errors
  X1  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24)), innov=e[,1])
  X2  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24)), innov=e[,2])

  #-- They seem to be correlated (because they are) --
  plot(X1, type="l", main=paste("cor=", round(cor(X1,X2), 4) ));  lines(X2, col="red")

  cor(X1,X2)
## [1] 0.6630245
  #-- Monte Carlo simulation shows --
  COR2 = 0
  for (i in 1:1000) {

    e = mvrnorm(n = 250, mu=c(0,0), Sigma=matrix(c(1,.7,.7,1), 2,2))
    X1  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24)), innov=e[,1])
    X2  = arima.sim(n = 250, list(ar = c(0.8), ma = c(-.24)), innov=e[,2])
    COR2[i] = cor(X1,X2)

  }
  hist(COR2, main="Cor of 1000 pairs of Dep St. TS", xlim=c(-1,1) )

When the two TS are stationary, there’s no problem using correlation as a guide for dependence.

Same goes for regressing two stationary TS.

3. Non-stationary TS

What happens when two non-stationary time series are regressed?

Simulate with two independent ARIMA(1,1,1) ?

Independent

    #-- X1 and X2 are both ARMA(1,1), but independent  --
    set.seed(63451)
    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))

    plot(X1, type="l", ylim=c(min(c(X1,X2)), max(c(X1,X2))),
        main=paste("cor=", round(cor(X1,X2), 4) ))
    lines(X2, col="red")

    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    plot(X1, type="l", ylim=c(min(c(X1,X2)), max(c(X1,X2))),  main=paste("cor=", round(cor(X1,X2), 4) ))
    lines(X2, col="red")

    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    plot(X1, type="l", ylim=c(min(c(X1,X2)), max(c(X1,X2))),  main=paste("cor=", round(cor(X1,X2), 4) ))
    lines(X2, col="red")

Repleat 1000 times and plot histogram of correlations

  #-- Repeat above in lopp to see overall behavior --
  n = 250
  COR4 = 0
  for (i in 1:1000) {
    X1  = arima.sim(n = n, list(order=c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = n, list(order=c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    COR4[i] = cor(X1,X2)
  }

  hist(COR4, main="Cor of 1000 pairs of Ind ARIMA(1,1,1)", xlim=c(-1,1) )

  quantile(COR4, c(.05, .95))
##         5%        95% 
## -0.7918580  0.8016304

This is the cause of what’s called Supurious Regression.

4. Spurious Regression

When you have two independent non-stationary time series, sample correlation tends to be high.

When correlation is high, regression parameters will be deemed significant.z

Can’t trust the result of regression, when two time series are non-stationary.

    #-- X1 and X2 are both ARMA(1,1), but independent  --
    set.seed(12342)
    X1  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))

    Reg1 = lm(X1~X2)
    summary(Reg1)
## 
## Call:
## lm(formula = X1 ~ X2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.726 -10.230   0.449   9.224  25.776 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -23.78038    1.51373  -15.71   <2e-16 ***
## X2            0.25396    0.02073   12.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13 on 249 degrees of freedom
## Multiple R-squared:  0.3761, Adjusted R-squared:  0.3736 
## F-statistic: 150.1 on 1 and 249 DF,  p-value: < 2.2e-16
    #-- They have nothing to do with each other, and correlation says so --
    plot(X1, type="l", ylim=c(min(c(X1,X2)), max(c(X1,X2))), main=paste("cor=", round(cor(X1,X2), 4) ))
    lines(X2, col="red")


I(\(d\)) notation

\(d\) is called Order of Integration

Series is called I(\(d\)) of taking difference \(d\) times make it stationary.

Notation: \[\begin{align*} X_t \sim I(1) &:= X_t \mbox{ is integrated series of order 1} \\\\ &= \hspace{3mm} \bigtriangledown X_t \mbox{ is stationary } \end{align*}\]

e.g. If \(X_t\) is ARIMA(2,1,2), then \(X_t \sim I(1)\).

In other words, ARIMA(2,1,2) has order of integration of 1.

Dependent

When two non-stationary series have actually something to do with each other, they may look like this:

    #- When two non-stationary series have relationship
    set.seed(45363)
    X11  = arima.sim(n = 250, list(order = c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X12  = 5 + .5*X11 + rnorm(251, 0, 5)

    plot(X11, type="l", ylim=c(min(c(X11,X12)), max(c(X11,X12))),
        main=paste("cor=", round(cor(X11,X12), 4) ))
    lines(X12, col="red")


Engle-Granger Method

Let’s check if the residual from regression is stationary or not.

  Reg11 = lm(X11~X12)
  summary(Reg11)
## 
## Call:
## lm(formula = X11 ~ X12)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.4646  -6.6099   0.1719   5.6251  28.6043 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.56677    0.76451  -16.44   <2e-16 ***
## X12           1.81535    0.03607   50.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.936 on 249 degrees of freedom
## Multiple R-squared:  0.9105, Adjusted R-squared:  0.9101 
## F-statistic:  2533 on 1 and 249 DF,  p-value: < 2.2e-16
  plot(X11, type='o')
  lines(Reg11$coef[1] + X12*Reg11$coef[2], col="red")

  plot(Reg11$residuals, type="o", main="Difference between Blk and Red")
  abline(h=0)

  source('https://nmimoto.github.io/R/TS-00.txt')
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
  St.test11 = Stationarity.tests(Reg11$resid)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Warning in adf.test(A): p-value smaller than printed p-value
## Warning in pp.test(A): p-value smaller than printed p-value
## Warning in kpss.test(A): p-value smaller than printed p-value
  St.test11
##        KPSS  ADF   PP
## p-val: 0.01 0.01 0.01
  cor(X11,X12)
## [1] 0.9541927


Testing E-G method on Independent Series

If we use E-G method on two independent ARIMA(1,1,1), how many can we catch?

  options(warn=-1)  # supress warinings for now

  COR1  = 0
  ST.T1 = 0
  for (i in 1:1000) {

    X1  = arima.sim(n = 250, list(order=c(1,1,1), ar = c(0.8), ma = c(-.24) ))
    X2  = arima.sim(n = 250, list(order=c(1,1,1), ar = c(0.8), ma = c(-.24) ))

    Reg1    = lm(X1~X2)
    COR1[i] = cor(X1,X2)
    ST.T1[i] = Stationarity.tests(Reg1$resid)[2]  #- use adf test p-val

  }

  hist(COR1, main=paste("mean = ", round(mean(COR1), 4) ), xlim=c(-1,1) )

  head( cbind(COR1, ST.T1) )
##             COR1 ST.T1
## [1,] -0.03924636 0.437
## [2,] -0.59196908 0.717
## [3,]  0.40226567 0.712
## [4,] -0.58573970 0.634
## [5,] -0.03836368 0.560
## [6,]  0.80863980 0.594
  coint = (ST.T1 < .05)   #- pick out index that had small ADF value
  head( cbind(COR1[coint], ST.T1[coint]) )
##            [,1]  [,2]
## [1,] -0.2231585 0.019
## [2,] -0.2273354 0.010
## [3,]  0.5921759 0.011
## [4,]  0.6621029 0.020
## [5,]  0.6537126 0.010
## [6,] -0.4000309 0.010
  sum(coint)   #- number of spurious reg not caught by E-G method (out of 1000)
## [1] 70
  options(warn=0)  # Turn warinings back on


Cointegration

There are many other economic theory that implies cointegrated relationship

Money Demand Model

Permanent Income Model

Unbiased Forward Rates Hypothesis

Fisher Equation


5. ARMAX model

ARMAX(2,1) is the model \[ Y_t = \mathbf \beta \mathbf X_t + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \epsilon_t + \theta_1 \epsilon_{t-1} \]

Don’t confuse this with additive model \[\begin{align*} K_t \hspace{3mm} &= \hspace{3mm} \mathbf \beta \mathbf X_t + Y_t \\ & Y_t \hspace{3mm} \sim \hspace{3mm} ARMA(2,1) \end{align*}\]

If we write ARMAX using backwards operator, \[\begin{align*} Y_t &= \mathbf \beta \mathbf x_t + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \epsilon_t + \theta_1 \epsilon_{t-1} \\ \\ \Phi(B) Y_t &= \mathbf \beta \mathbf x_t + \Theta(B) \epsilon_t \\ \\ Y_t &= \frac{ \mathbf \beta \mathbf x_t }{ \Phi(B) } + \frac{ \Theta(B) }{ \Phi(B) } \epsilon_t \end{align*}\] which is very hard to interpret.

Currently, there is no package that directly deal with ARMAX model.