Además de realizar inferencias sobre las variables que causan a otra, el modelo de regresión se puede emplear para realizar predicciones.

library(forecast, quietly=TRUE)
library(ggplot2, quietly=TRUE)
library(dplyr, quietly=TRUE)
library(tseries, quietly=TRUE)
data(economics)
datos <- economics[,2:3]
datos_stationary = datos %>%
  lapply(log) %>%
  lapply(diff) %>% 
  as.data.frame() %>%
  ts(start=c(start=c(1967,7), end=c(2015, 03)), deltat = 1/12)#Log-difenciar los datos y hacerlos series temporales
autoplot(datos_stationary) + theme_classic()

lapply(datos_stationary, kpss.test) #Aplicar test de estacionareidad

$pce

    KPSS Test for Level Stationarity

data:  X[[i]]
KPSS Level = 4.4404, Truncation lag parameter = 6, p-value = 0.01


$pop

    KPSS Test for Level Stationarity

data:  X[[i]]
KPSS Level = 1.5535, Truncation lag parameter = 6, p-value = 0.01

datos_st <- datos_stationary %>%
  lapply(diff) %>%
  as.data.frame() %>%
  ts(start=c(start=c(1967,7), end=c(2015, 02)), deltat = 1/12)  #Volver a log-diferenciar
autoplot(datos_st) + theme_classic()

lapply(datos_st, kpss.test) #Test

$pce

    KPSS Test for Level Stationarity

data:  X[[i]]
KPSS Level = 0.0063759, Truncation lag parameter = 6, p-value = 0.1


$pop

    KPSS Test for Level Stationarity

data:  X[[i]]
KPSS Level = 0.012085, Truncation lag parameter = 6, p-value = 0.1

modelots <- tslm(pce ~ pop, data=datos_st) #Modelo
nuevosdatos = data.frame(pop=c(1.406604e-05, -2.109815e-05, 1.196598e-04, -2.567329e-05)) #Generar vector de variable independiente para predecir.
pred <- forecast(modelots, newdata=nuevosdatos) #Predecir
autoplot(pred) + theme_classic()

Aunque, conforme se llevan a cabo transformaciones de los datos, tanto las predicciones como los coeficientes son cada vez menos interpretables.

Dummies estacionales

Para evitar tener que hacer diferenciación estacional, se pueden aplicar variables ficticias o dummies.

$Y_{t} = β_{0} + \sum_{i = 1}^{k} β_{i} X_{i, t} + \sum_{j = 1}^{s} α_{j} d_{j, t} + ϵ_{t}$

donde cada variable $d$ es una variable que toma valor 1 si es el periodo del año al que corresponde y 0 si no. La clave está en añadir una variable para cada estación del año $1, \dots, s$ menos una, para evitar una combinación lineal perfecta de los predictores (trampa de las variables ficticias). Los coeficientes de las variables recogerán el efecto del componente estacional. Si a la serie original, se le sustrae este efecto, se tiene lo que se conoce como una serie desestacionalizada.

Si una serie posee una tendencia y un componente estacional, pero no un componente estocástico, es posible hacerla estacionaria mediante la inclusión de dummies estacionales, y una función de la tendencia.

data("gas")
modelogas <- tslm(gas ~ trend + season, data=gas) #Modelo con tendencia determinista y dummies estacionales
pred.gas <- forecast(modelogas)
autoplot(pred.gas) + theme_classic()

Predicción de escenarios

library(Ecdat, quietly=TRUE)
data("Klein")
Kleints <- as.ts(Klein, start=1920, end=1940)
Kleints <- lapply(Kleints, diff) %>% as.data.frame() %>% as.ts(start=1920, end=1939)
autoplot(Kleints[,c(8,9)]) + theme_classic()

kpss.test(diff(Kleints[,8]))#Test de estacionareidad


    KPSS Test for Level Stationarity

data:  diff(Kleints[, 8])
KPSS Level = 0.42086, Truncation lag parameter = 2, p-value = 0.06817

kpss.test(diff(Kleints[,9])) #Test de estacionareidad


    KPSS Test for Level Stationarity

data:  diff(Kleints[, 9])
KPSS Level = 0.14562, Truncation lag parameter = 2, p-value = 0.1

gobmod <- tslm(govspend ~ taxe, data=Kleints) #Modelo de series temporales
gobmod


Call:
tslm(formula = govspend ~ taxe, data = Kleints)

Coefficients:
(Intercept)         taxe  
     0.4336       0.2798

h = 2 #Establecer un horizonte temporal
Kleints #Visualizar valores razonables para los impuestos

Time Series:
Start = 1 
End = 21 
Frequency = 1 
   cons profit privwage  inv lcap  gnp pubwage govspend taxe
 1  2.1   -0.3     -3.3 -2.9  2.7  0.7     0.5      1.5  4.3
 2  3.1    4.5      3.8  2.1 -0.2  4.5     0.2     -0.7 -3.8
 3  4.2    1.5      4.8  3.3  1.9  7.1     0.0     -0.4  0.8
 4  1.4    1.0     -0.2 -2.2  5.2 -0.1     0.2      0.7 -0.9
 5  2.0    0.7      1.5  2.1  3.0  3.9     0.1     -0.2  1.7
 6  2.5   -0.5      2.0  0.5  5.1  3.0     0.1      0.0  1.5
 7  1.1    0.2      0.5 -1.4  5.6  0.4     0.3      0.7 -0.3
 8  1.1    1.3      1.3 -1.2  4.2  0.1     0.1      0.2 -2.5
 9  0.5    0.6      2.1  2.1  3.0  2.5     0.3     -0.1 -0.2
10 -2.8   -6.1     -3.4 -4.1  5.1 -5.8     0.2      1.1  3.7
11 -4.1   -4.2     -3.4 -4.4  1.0 -7.8     0.6      0.7 -0.2
12 -5.3   -4.4     -5.5 -2.8 -3.4 -9.1     0.5     -1.0  0.8
13  0.9    4.2     -0.5  1.1 -6.2  0.8     0.3     -1.2 -2.9
14  2.2    1.1      2.1  2.1 -5.1  4.6     0.4      0.3  1.4
15  2.6    1.7      2.6  1.7 -3.0  4.7     0.1      0.4  0.4
16  6.4    3.6      3.6  3.4 -1.3  8.3     1.3     -1.5  1.1
17  1.0   -0.3      4.2 -0.1  2.1  2.3    -0.7      1.4 -1.6
18 -1.2   -2.0     -2.8 -3.9  2.0 -4.1     1.0      1.0  0.7
19  4.1    3.7      3.4  3.2 -1.9  8.6     0.1      1.3  1.5
20  3.4    2.1      3.4  2.0  1.3  6.2     0.2      0.8  0.7
21  4.7    2.4      8.3  1.6  3.3 12.7     0.5      6.4  2.0

valores.opt <- data.frame(taxe=c(10, 12)) #Escenario optimista
valores.pes <- data.frame(taxe=c(-6, -7)) #Escenario pesimista
pred.opt <- forecast::forecast(gobmod, newdata=valores.opt) #Predicción optimista
pred.pes <- forecast::forecast(gobmod, newdata=valores.pes) #Predicción pesimista
autoplot(Kleints[,8]) + xlab("") + ylab("") + autolayer(pred.opt, PI=T, series="Optimista") +  autolayer(pred.pes, PI=T, series = "Pesimista") + guides(colour=guide_legend("Escenario")) + theme_classic()

Ejercicio: Con solo una línea de código, prediga el comportamiento de la serie AirPassengers a través de su tendencia y dummies estacionales. ¿Es una predicción razonable?