Workshop 2, Time Series

Abstract: In this workshop we introduce the concepts of stationarity, unitary root and cointegration of time-series variables.

##Stationarity of time series variables

In this section we will work with historical sales of consumer products to examine the concept of stationarity in time-series.

Introduction

In our previous workshop we learned about the random walk model applied to finance. We learned that the logarithm of the stock prices behaves similar to a random walk model with a positive drift (ϕ0>0).

As we learned, a random walk model can be expressed in the following equation:

\[Y_{t}=\varphi _{0}+Y_{T-1}+\varepsilon _{t}\] Following this model we did a simulation for the logarithm of the S&P500, and we found that most of the time the simulations of a random walk model behaves similar to the real log of the S&P500 index.

We learned that if ϕ0>0, then the series will have a growing trend over time. We also learned that the higher the standard deviation of the daily shock, the more likely that the series will radically move up and down. The standard deviation of the shock is actually the daily volatility of the daily returns of the series.

An important insight about the random walk model is that the random walk model is a non-stationary series. What does this mean?

In short, a stationary series is a series that its mean in any time period is about the same (no matter which time periods we select), and its standard deviation is about the same for any time period. There is another condition of stationarity related to auto-covariance of the series, but we will leave the concept of auto-covariance for a later workshop.

Then, the random-walk model will always behave as a non-stationary series.

Let’s analyze what part of the random walk equation makes the series stationary. We can re-write the previous random walk equation as follows:

\[Y_{t}=\varphi _{0}+\phi _{1}*Y_{T-1}+\varepsilon _{t}\] If the coefficient ϕ1=1, then this equation becomes the random walk equation. When ϕ1=1 the Yt series becomes non-stationary.

If ϕ1<1, then Yt becomes a stationary series.

Let’s work with an example.

##Data collection and data management

We will work with an online dataset that contains sales data of different consumer products. The data was created using real products, but the names of the products were changed due to privacy issues.

Download the excel file from a web site:

download.file("http://www.apradie.com/ec2004/salesbrands3.xlsx", "salesbrands3.xlsx", mode="wb")

This dataset has sales information of several brands of consumer products that belong to one category (for example, the category can be Cereals or Candies). The content of this dataset is very similar to a real-world category. Consumer products usually show some seasonality and trend over time. We will learn how to handle these features later.

Also, you have to install the readxl package. This package reads data from Excel files. We load this package:

library(readxl)

Now, we import the Excel file into an R dataframe:

sales_brands <- read_excel("salesbrands3.xlsx", sheet = "Resultado")

# Always familiarize with the data
head(sales_brands, 5)

As you can see, we have imported a monthly dataset where the variable date shows the first day of the month in a yyyy-mm-dd format.

In order to perform some time-series analysis, it is recommended to use xts objects. We will construct a xts object using the xts() method. The objects of xts class are composed by 2 objects: the core data and the index (class Date, POSIX, etc).

# Load library xts. Install the package if you haven't already.
library(xts)

# We will start by turning sales_brands into a data frame 
# because data frames are easier to manipulate
sales_brands.df <- as.data.frame(sales_brands)

# Divide the index and the core data
fecha <- sales_brands.df$date
sales_brands.df <- sales_brands.df[,-1]

# Join both objects using xts() and assign the resulting object to sales.xts
sales.xts <- xts(x = sales_brands.df, order.by = fecha)

head(sales.xts, 5)

timezone of object (UTC) is different than current timezone ().

           qbrand1tons salesbrand1 qbrand2tons salesbrand2 qbrand3tons salesbrand3 qbrand4tons salesbrand4 qbrand5tons salesbrand5
2015-01-01      106.62    60780.72      195.46    60700.51      293.10    75196.35      185.91    48228.51       28.94     2295.50
2015-02-01      105.82    60411.47      216.09    68289.15      494.46   105010.90      149.30    39922.14       24.66     1973.43
2015-03-01       88.82    53174.89      220.15    62529.54      356.88    90902.84      144.65    39788.22       23.75     1903.24
2015-04-01       86.69    50861.69      170.74    50575.68      291.07    72853.34      161.43    42692.59       21.00     1684.28
2015-05-01       87.23    51014.23      177.17    54327.48      318.09    75277.10      141.09    37167.52       23.44     1883.41
           qbrand6tons salesbrand6 tonstotal salestotal
2015-01-01      338.38    55561.26   1148.41   302762.8
2015-02-01      284.52    51469.74   1274.85   327076.8
2015-03-01      251.34    50395.48   1085.59   298694.2
2015-04-01      230.29    49124.61    961.22   267792.2
2015-05-01      250.45    46906.63    997.47   266576.4

# Do you notice any difference between data frames and xts objects?

#delete variable fecha
rm(fecha)

##Stationary vs Non-stationary

We can have an overview of the behavior of sales of these consumer products using the plot command. Let’s focus on brand 2 and 6. To see volume sales (the number of tons sold) of brand 2:

plot(sales.xts$qbrand2tons)

It seems that this product has a growing trend over time, except for the last year.

We can graph consumer brand 6:

plot(sales.xts$qbrand6tons)

This brand has a clear growing trend over time. As we can see, many consumer products have a growing or a declining trend. And other products might seem to be stationary over time, with sales over an average over time.

A series with a growing or declining trend is non-stationary. Most of the variables that represent sales of any product can be treated as non-stationary. Then, a rule of thumb when we work with business variables such as volume sales, value sales, cost of good sold, etc is to transform the series into a stationary series.

A SERIES WITH A GROWING OR DECLINING TREND IS NS.

When you work with sales data, I strongly recommend to create the first difference of the logarithm of the variable. This is actually its percentage change (in continuously compounded). If you do this process to most of real-world sales variable the result will be stationary most of the time (about 95% of the variables). The other reason why it is recommended to do this transformation is that the first difference of the log of a variable is actually the % change of the variable from one period to the next (in continuously compounded percentage). Then, we can do this transformation as follows:

THE FIRST DIFFERENCE OF THE LOG OF A VARIABLE IS ACTUALLY THE % CHANGE OF THE VARIABLE FROM ONE PERIOD TO THE NEXT.

We first create the natural log of each sales variable (in tons). We can apply the log function to the xts object to get the log of all series:

ln_sales <- log(sales.xts)

The new object has the same structure as sales.xts but the columns contain the natural logarithm of the data. This is due to R vectorization, a key feature of R.

Now, I can work with the first difference month-by-month using the diff() function for each variable. For example, to graph the first difference of the log of volume sales for both brands :

plot(diff(ln_sales$qbrand2tons), type = "l", col = "red")

lines(diff(ln_sales$qbrand6tons), type = "l", col = "blue")

The first difference of the log value in R is diff(ln_sales$qbrand2tons). The diff() function generates on the fly this first difference as:

diff(ln_salesqbrand2tons∗) = ∗lnsalesqbrand2tons - lag(ln_sales$qbrand2tons, 1)

Remember that the dollar sign operator ($) is used to access/refer to the columns of a data frame or xts object.

We can see now that the 2 series look like stationary since their means tend to be in the same level.

R performs the first difference calculation temporarily to do a graph, and then, this variable is not saved in the dataset or environment. In other words, we do not need to generate the first difference as another variable (column), but we can if we need to.

Another transformation with monthly data is **annual percentage change*, which can be calculated as the difference between the log value of the series and its log value 12 months ago. We can check this difference using the diff() function with the lag argument set to 12:

plot(diff(ln_sales$qbrand2tons, lag = 12), type = "l", col = "red")

lines(diff(ln_sales$qbrand6tons, lag = 12), type = "l", col = "blue")

In this case, it is difficult to evaluate graphically whether both series are stationary. We need to do a statistical test to check whether a series can be treated as stationary or not. Some series might not look stationary for our eyes, but we need to run statistical tests to decide whether to treat the series stationary or not.

There is a statistical way to evaluate whether these series are stationary. We need a statistical confirmation in addition to the visual tools. The Dicky-Fuller test is one of the most used statistical tests for stationary. This test is called a unit root test. The Dicky-Fuller test assumes that the series is non stationary. In other words, it assumes that the series follows a random walk with a unit root (ϕ1=1).

DICKY-FULLER ASSUMES THAT THE SERIES IS NON STATIONARY.

You need to install the package tseries to do the following analysis.

We can run the Dicky-Fuller test for the seasonal difference of the product 6 log sales as follows:

library(tseries)
adf.test(na.omit(diff(ln_sales$qbrand2tons, lag = 12)), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(diff(ln_sales$qbrand2tons, lag = 12))
Dickey-Fuller = -3.574, Lag order = 3, p-value = 0.04701
alternative hypothesis: stationary

The adf.test() function from the tseries package will do a Dickey-Fuller test if we set lags (k) equal to 0. Use the default value of k.

The null hypothesis of the Dicky-Fuller test is that the series is non-stationary. Actually, the *null-hypothesis of this test is that the ϕ1=1, which implies that the series is non-stationary. This is the reason why the Dicky-Fuller test is called a unit-root** test.

HO = ϕ1=1 -> NON STATIONARY

As any statistical test, if the p-value of the test is less than 0.05, then we have statistical evidence (at the 95% confidence level) to REJECT the null-hypothesis.

This is another way to do the same Dickey Fuller test:

# I create a variable for the seasonal difference of the log sales:
s12.lnbrand2tons <- na.omit(diff(ln_sales$qbrand2tons, lag = 12))
# I run Dicky-Fuller test for this seasonal difference:
adf.test(s12.lnbrand2tons, alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  s12.lnbrand2tons
Dickey-Fuller = -3.574, Lag order = 3, p-value = 0.04701
alternative hypothesis: stationary

The Dicky Fuller test reports a p-value less than 0.05, indicating that we can reject the null hypothesis that states a unit root for these series (ϕ1=1). In other words, since the p-value is less than 0.05 we can say that there is statistical evidence to conclude that the seasonal difference of the log series (the annual % change) (s12.lnbrand2tons) is stationary.

Then for the variable s12.lnbrand2tons we can say that there is statistical evidence to say that the series is stationary.

RUN A DICKY-FULLER TEST FOR THE ANNUAL % DIFFERENCE OF BRAND 6, AND PROVIDE AN INTERPRETATION

s12.lnbrand6tons <- na.omit(diff(ln_sales$qbrand6tons, lag = 12))
adf.test(s12.lnbrand6tons, alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  s12.lnbrand6tons
Dickey-Fuller = -2.3423, Lag order = 3, p-value = 0.4387
alternative hypothesis: stationary

FOR THE VERY 6TH BRAND, WE MAY INTERPRET THAT THERE IS STATISTICAL EVIDENCE TO CONCLUDE THAT THE SEASONAL DIFFERENCE OF THE ANNUAL % CHANGE (s12.lnbrand6tons) IS NON STATIONARY. IN THIS CASE, THE DICKY FULLER TEST REPORTS A P-VALUE HIGHER THAN 0.05, INDICATING US THAT WE CAN ACCEPT THE NULL HYPOTHESIS THAT STATES HO = ϕ1=1 -> NON STATIONARY.

For the case of monthly times series variables like sales of products or sales of a company, I strongly recommend to work with the seasonal difference of the log of sales (s12.lnbrand2tons), which is the annual % change (in continuously compounding) month by month. The reason is because most of consumer products have a seasonal pattern. We can do this ONLY if the seasonal difference of the log is STATIONARY.

For these 2 products, then, I can use the seasonal difference of the log of sales to design a time-series model.

##Spurious regression

When we want to examine the relationship between two non-stationary variables by running a regression model, we have the risk to end up with a non-valid - spurious - regression. Before we understand why a regression model can be spurious, we start with and example using 2 real-world variables.

Install the wbstats package using RStudio. This package was written by the World Bank, and it has functions to download data of all countries around the world. The function wb downloads hundreds of time-series variables that the World Bank tracks for all countries. If you want to know more about this package, you can check its documentation in the cran web site (https://cran.r-project.org/web/packages/wbstats/wbstats.pdf)

We will download the infant mortality rate and the exports for Mexico. It is supposed that these variables have nothing in common, so we would not expect a significant relationship.

library(wbstats)

# Mexico - Infant mortality
infantm<-wb_data(indicator = c("SP.DYN.IMRT.IN"), 
      country="MEX", start_date = 1980, end_date = 2020)

# Mexico - Export value
exports<-wb_data(indicator = c("TX.VAL.MRCH.XD.WD"), 
      country="MEX", start_date = 1980, end_date = 2020)

The wb function brings a data frame with the requested data. We can plot the data to have an idea of these 2 variables:

plot.ts(infantm$SP.DYN.IMRT.IN)

plot.ts(exports$TX.VAL.MRCH.XD.WD)

Now run a regression using these series. Report the result of the regression. Did you find significant relationship? is your result what you expected?

COMO PODEMOS VER, CLARAMENTE EXISTE UNA RELACIÓN INVERSA ENTRE LAS VARIABLES; MIENTRAS QUE UNA CRECE, LA OTRA DECRECE. SIN EMBARGO, EL HECHO DE QUE ESTAS VARIABLES ESTÉN RELACIONADAS, NO NECESARIAMENTE NOS INDICA QUE ESTÁN COEINTEGRADAS, Y PARA ELLO EMPLEAREMOS EL DICKEY FULLER TEST.

m1 <- lm(exports$TX.VAL.MRCH.XD.WD ~ infantm$SP.DYN.IMRT.IN)
summary(m1)

Call:
lm(formula = exports$TX.VAL.MRCH.XD.WD ~ infantm$SP.DYN.IMRT.IN)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.331 -39.723  -7.456  37.848  77.194 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            275.1783    16.0408   17.16  < 2e-16 ***
infantm$SP.DYN.IMRT.IN  -6.1847     0.5327  -11.61  4.6e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 41.65 on 38 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.7801,    Adjusted R-squared:  0.7743 
F-statistic: 134.8 on 1 and 38 DF,  p-value: 4.602e-14

RESEARCH ABOUT SPURIOUS REGRESSION. With your words briefly explain in which cases you can end up with a spurious regression

UNA RELACIÓN ESPURIA HACE REFERENCIA A LA APARIENCIA QUE EXISTE EN UNA RELACIÓN DE CAUSALIDAD ENTRE VARIABLES CUANDO EN REALIDAD ESTA NO EXISTE. ES DECIR, QUE AUNQUE PARECIERA QUE ALGUNAS VARIABLES ESTÁN CORRELACIONADAS, REALMENTE NO HAY UNA RELACIÓN CAUSAL ENTRE ESTAS, POR LO TANTO ES ESPURIA. POR EJEMPLO, LAS VENTAS DE NIEVES EN LA CIUDAD DE QUERÉTARO Y LAS MUERTES POR COVID-19 EN ESTADOS UNIDOS. AUNQUE EXISTIERA LA POSIBILIDAD DE DETERMINAR CIERTA RELACIÓN, EN CASO DE QUE EXISTA, SERÍA ESPURIA, PORQUE REALMENTE NO TIENEN NADA QUE LA UNA CON LA OTRA.

##Cointegration

When two non-stationary series are strongly related in the long run, then we say that both series are cointegrated. The question is how strong this long-term relationship needs to be in order to consider 2 non-stationary series as cointegrated.

To do a cointegration statistical test of 2 non-stationary series, we need to do the following:

1. Run a linear regression of the 2 non-stationary series

2. Get the residuals of this regression

3. Do a Dicky-Fuller test to check whether this regression residuals (errors) is a stationary series. If the p-value<0.05, then we can say that there is statistical evidence to consider both non-stationary series as cointegrated.

If 2 non-stationary series are NOT cointegrated, then its linear regression becomes a spurious regression. In other words, when 2 non-stationary series are not cointegrated, the result of a linear regression using these 2 variables will not be reliable or valid.

Are the infant mortality and export series cointegrated? Is the regression spurious or valid? Run and interpret the corresponding test

Here is a code to help you respond this question.

library(tseries)
# The errors/residual of the regression are stored in the residuals
#   attribute of the m1 model:
adf.test(m1$residuals, k = 0)

    Augmented Dickey-Fuller Test

data:  m1$residuals
Dickey-Fuller = -1.5881, Lag order = 0, p-value = 0.7351
alternative hypothesis: stationary

Tras haber empleado el Dickey-Fuller test, podemos determinar que la relación entre la mortalidad infantil y las exportaciones es espuria. Ya que, aunque acertadamente podemos apreciar una relación inversa; no hay suficiente soporte estadístico para determinar que ambas variables estén correlacionadas. Simplemente es una relación espuria. Con el valor p, en este caso de 73.51%, terminamos por concluir que la Ho es correcta, mientras que la Ha es rechazada. Es decir, si nos fiamos de dicha relación para emplear algun estudio, existe un 73.51% de equivocación en nuestros hallazgos.

##Cointegration between Financial series

Using daily data of Mexico IPCyC market index and the S&P 500, examine whether two series are cointegrated. Generate an index for each instrument. To do these indexes, create a variable that represents how 1.00 peso or 1.00 dollar invested in each instrument would be changing over time.

1. From Jan 1, 2015 to Oct 2, 2017.
2. From Oct 3, 2017 to Feb 28, 2018

For both cases, run a cointegration test and INTERPRET your results.

rm(list=ls())
options(scipen=999)

library(quantmod)

getSymbols(Symbols= c("^MXX", "^GSPC"), from="2015-01-01", to= "2017-10-02", 
           periodicity = "daily", src = "yahoo")

[1] "^MXX"  "^GSPC"

prices = Ad(merge(MXX,GSPC))
prices = na.omit(prices)

MXX1 = as.numeric(MXX$MXX.Adjusted[1])
GSPC1 = as.numeric(GSPC$GSPC.Adjusted[1])

ONEMXP = prices$MXX.Adjusted / MXX1
ONEUSD = prices$GSPC.Adjusted / GSPC1

plot(ONEMXP$MXX.Adjusted)

plot(ONEUSD$GSPC.Adjusted)

m3 <- lm(ONEMXP$MXX.Adjusted ~ ONEUSD$GSPC.Adjusted)
summary(m3)

Call:
lm(formula = ONEMXP$MXX.Adjusted ~ ONEUSD$GSPC.Adjusted)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.06656 -0.01835  0.00244  0.01944  0.06810 

Coefficients:
                     Estimate Std. Error t value            Pr(>|t|)    
(Intercept)           0.35668    0.01383   25.79 <0.0000000000000002 ***
ONEUSD$GSPC.Adjusted  0.69852    0.01311   53.29 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02702 on 672 degrees of freedom
Multiple R-squared:  0.8086,    Adjusted R-squared:  0.8083 
F-statistic:  2840 on 1 and 672 DF,  p-value: < 0.00000000000000022

library(tseries)
adf.test(m3$residuals, k = 0)

    Augmented Dickey-Fuller Test

data:  m3$residuals
Dickey-Fuller = -3.5752, Lag order = 0, p-value = 0.03508
alternative hypothesis: stationary

ES BIEN SABIDO QUE LA ECONOMÍA MEXICANA DEPENDE DE LA ECONOMÍA ESTADOUNIDENSE EN GRAN PORCENTAJE; NO OBSTANTE, CON EL FIN DE DETERMINAR QUE DICHA CORRELACIÓN NO SEA ESPURIA, DECIDIMOS CORRER EL DICKEY-FULLER TEST Y LOGRAMOS DETERMINAR QUE EFECTIVAMENTE AMBOS ÍNDICES ESTÁN COINTEGRADOS, ASÍ COMO LA EXISTENTE CORRELACIÓN ENTRE ESTAS VARIABLES. CON UN VALOR P DE 3.5%, TAMBIÉN COMPROBAMOS ESTADÍSTICAMENTE QUE AMBAS SERIES SON NO ESTACIONARIAS Y UNA VEZ MÁS, QUE ESTÁN COINTEGRADAS LAS VARIABLES.

rm(list=ls())
options(scipen=999)

library(quantmod)

getSymbols(Symbols= c("^MXX", "^GSPC"), from="2017-10-03", to= "2018-02-28", 
           periodicity = "daily", src = "yahoo")

[1] "^MXX"  "^GSPC"

prices2 = Ad(merge(MXX,GSPC))
prices2 = na.omit(prices2)

MXX2 = as.numeric(MXX$MXX.Adjusted[1])
GSPC2 = as.numeric(GSPC$GSPC.Adjusted[1])

ONEMXP2 = prices2$MXX.Adjusted / MXX2
ONEUSD2 = prices2$GSPC.Adjusted / GSPC2

plot(ONEMXP2$MXX.Adjusted)

plot(ONEUSD2$GSPC.Adjusted)

m4 <- lm(ONEMXP2$MXX.Adjusted ~ ONEUSD2$GSPC.Adjusted)
summary(m4)

Call:
lm(formula = ONEMXP2$MXX.Adjusted ~ ONEUSD2$GSPC.Adjusted)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.037454 -0.014008 -0.002567  0.017749  0.040777 

Coefficients:
                      Estimate Std. Error t value             Pr(>|t|)    
(Intercept)            0.79716    0.05928  13.447 < 0.0000000000000002 ***
ONEUSD2$GSPC.Adjusted  0.16206    0.05640   2.874              0.00501 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02006 on 95 degrees of freedom
Multiple R-squared:  0.07997,   Adjusted R-squared:  0.07028 
F-statistic: 8.257 on 1 and 95 DF,  p-value: 0.005008

library(tseries)
adf.test(m4$residuals, k = 0)

    Augmented Dickey-Fuller Test

data:  m4$residuals
Dickey-Fuller = -2.0418, Lag order = 0, p-value = 0.5592
alternative hypothesis: stationary

EN ESTE SEGUNDO PERIODO DE TIEMPO, CURIOSAMENTE OBTENEMOS DISTINTOS RESULTADOS. LA REGRESIÓN LINEAL ESTÁ CLARA, SI EXISTE UNA RELACIÓN POSITIVA ENTRE AMBAS VARIABLES Y EL VALOR P (DE LA REGRESIÓN LINEAL) NOS RECTIFICA CON CERTEZA SU ALTA RELACIÓN. NO OBSTANTE, TRAS HABER CORRIDO EL DICKEY FULLER TEST, PODEMOS OBSERVAR CÓMO ESTE VALOR SUPERA POR MUCHO EL 5%; INDICANDONOS QUE NO EXISTE COINTEGRACIÓN ENTRE LAS VARIABLES Y QUE AMBAS SERIES SON ESTACIONARIAS.