Econometría II: Tendencia

Departamento de Economía

Carlos A. Yanes G.

2024-08-01

Paquetes con que se trabaja la sesión

R libraries

Los paquetes que vamos a utilizar y que se deben instalar para ser usados en clase son:

install.packages(c("pacman", "fpp2", "dynlm", "huxtable"))

Ejecución

No olvide que debe desde luego tener presente en su documento de R Markdown cargar estos paquetes.

library(pacman)
p_load(fpp2, dynlm, huxtable)

Datos

Cartera comercial de bancos

Paquete fpp2

  • Cuando usted carga library(fpp2) de (Hyndman and Athanasopoulos 2018). Este contiene elementos y/o funciones de otros paquetes como:
    • forecast Para pronostico simple
    • ggplot2 Gráficos como si hiciera uso de
    • Ademas de fma y expsmooth para muchas series

Modelos de Tendencia central

  • Estimamos a partir de la tendencia un pequeño pronostico de la serie de la cartera comercial.
    • Modelo de tendencia lineal \(Y_t=f(T)\)
    • Modelo de tendencia de polinomio superior
    • Modelo exponencial
    • Modelo Box Cox

Caso de estudio

Caso 2

Se solicita lo siguiente:

  1. Caracterización general de la base de datos.
  2. Trabajar con la tendencia y escoger la mejor especificación
  3. Hacer una predicción.

Solución en

Caso 2

Code 
# Cargar base de datos
library(readxl)
bsd2 <- read_excel("cartera.xlsx")
head(bsd2)
# A tibble: 6 × 2
  Fecha               Cartera
  <dttm>                <dbl>
1 2016-01-01 00:00:00 367822.
2 2016-02-01 00:00:00 372018.
3 2016-03-01 00:00:00 371504.
4 2016-04-01 00:00:00 373540.
5 2016-05-01 00:00:00 378923.
6 2016-06-01 00:00:00 381253.
Code 
# Cargar base de datos
library(readxl)
library(tidyverse)
library(fpp2)
datos <- read_excel("cartera.xlsx")
datos <- datos$Cartera
bsd2 <- ts(datos, start = c(2016, 1), frequency = 12)
bsd2
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2016 367822.1 372017.9 371504.5 373540.2 378922.9 381252.9 384245.6 385870.2
2017 389852.2 393366.3 394979.5 397288.4 399964.5 403662.3 403451.3 404165.9
2018 410308.2 412205.1 414766.3 416397.9 419080.9 420973.1 420294.1 421235.2
2019 433608.6 438173.4 441902.2 444898.5 449547.2 451756.4 454250.3 458674.4
2020 472036.3 478144.4 496248.9 499381.2 500477.1 501224.8 497737.8 492528.9
2021 485300.9 489968.3 492815.4 496422.9 500597.8 505475.2 510043.2 515702.4
2022 544147.8 554838.4 563242.2 573378.8 582534.2 593664.9 601094.8 608790.9
2023 633633.3 637485.1 636143.1 638791.7 640048.1 644309.1 644511.0 646039.0
2024 646722.0                                                               
          Sep      Oct      Nov      Dec
2016 387285.4 388795.6 392844.9 393753.7
2017 407058.0 408665.0 411022.7 413309.0
2018 423279.1 428474.8 434026.8 437514.9
2019 464159.0 466107.9 472113.6 472288.1
2020 490861.8 489415.5 488959.8 487225.6
2021 522266.9 527621.2 537485.1 543515.2
2022 616944.7 624696.6 630909.7 636903.2
2023 648621.0 648701.0 650135.0 649298.0
2024

  • El código para gráficar:
autoplot(bsd2) + ylab("$ miles de Millones") + xlab("Year") +
         ggtitle("Cartera Comercial de Bancos Colombianos 2016-2023")

  • autoplot Es una función o comando del paquete fpp2. Contiene elementos de ggplot2 quien es el mas importante a la hora de hacer gráficos en todo el ecosistema de

  • Lo importante es mirar siempre la tendencia de la serie.

Estacionalidad

Estacionalidad

ggseasonplot(bsd2, year.labels=TRUE, year.labels.left=TRUE) +
  ylab("$ Miles de millones") +
  ggtitle("Gráfico Estacional: Cartera comercial")
  • Si nota, automáticamente la parte de ggseaasonplot descompone la serie por periodo o fecha que lo tiene. Recuerde que para hacer eso el objeto que este trabajando debe necesariamente tener el index de ts.

Algunas funciones que trae consigo autoplot son:

  • seasonal() Para extraer el componente estacional
  • trendcycle() Para la parte de tendencia-ciclo
  • remainder() Extrae la parte del componente irregular.
  • seasadj() Ajusta la serie por sesión.

Estacionalidad

Estacionalidad

sea <- stl(bsd2, s.window="per")
autoplot(bsd2, series="Data") +
  autolayer(seasadj(sea), series="Ajuste Estacional")
  • La opción de s.windows puede entenderse de forma muy especial si se sigue el documento de (Cleveland et al. 1990). Las recomendaciones en ventanas periodicas como s.window="per" o como números como s.window=7.
  • Tenemos una serie que puede ser aditiva:

\[y_t= S_t+T_t+e_t\]

  • Si deseamos eliminar su estacionalidad es “fácil” decir que:

\[y_t - S_t = T_t + e_t\]

  • En caso tal la miraramos por el lado multiplicativo

\[\frac{y_t}{S_t} = T_t \times e_t\] Ya que:

\[y_t= S_t\times T_t\times e_t\]

Ciclo

Ciclo y tendencia

Code 
autoplot(bsd2, series="Data") +
  autolayer(trendcycle(sea), series="Tendencia") +
  autolayer(seasadj(sea), series="Ajuste Estacional") +
  xlab("Fecha") + ylab("$ Miles de Millones") +
  ggtitle("Cartera de Bancos comerciales") +
  scale_colour_manual(values=c("gray","blue","red"),
                     breaks=c("Data","Ajuste Estacional","Tendencia"))

Tendencia

Tendencia

  • Vamos a correr varios modelos con las tendencias a ver que ajuste tenemos. Será nuestro intento de primer pronostico

  • Podemos usar cualquier paquete de estimación. Sin embargo en esta parte para el resumen haremos uso entonces del paquete huxtable.

Modelo Tendencia lineal
Variables Parámetro SE t P-value
(Intercept) 337984.32 4448.849 75.97 0
trend(bsd2) 36970.09 945.963 39.08 0
Modelo Tendencia Cúbica
Variables Parámetro SE t P-value
(Intercept) 373580.083 5366.119 69.62 0.000
trend(bsd2) 14398.756 5661.322 2.54 0.013
I(trend(bsd2)^2) 1537.276 1606.378 0.96 0.341
I(trend(bsd2)^3) 167.452 129.337 1.29 0.199 
library(dynlm)
library(huxtable)
datos <- read_excel("cartera.xlsx")
bsd2 <- ts(datos, start = c(2016, 1), frequency = 12)
# Ajustar el modelo de regresión lineal solo con tendencia
m1 <- dynlm(Cartera ~ trend(bsd2), data = bsd2)
m2 <- dynlm(Cartera ~ trend(bsd2)+I(trend(bsd2)^2), data = bsd2)
m3 <- dynlm(Cartera ~ trend(bsd2)+I(trend(bsd2)^2)+I(trend(bsd2)^3), data = bsd2)
m4 <- dynlm(Cartera ~ log(trend(bsd2)), data = bsd2)
m5 <- dynlm(Cartera ~ I(1/trend(bsd2)), data = bsd2)
huxreg(m1,m2,m3,m4,m5, statistics=c("Observaciones" = "nobs", "R2 Ajustado" = "adj.r.squared"))

Primer pronostico 😲

Primer pronostico

Imaginemos un modelo de tendencia polinomica

  • Se toma la ecuación del modelo (3) y se reemplaza en el último valor o periodo tendencial. P.e:

\[\begin{aligned} Cartera_t=& \; \alpha_0+\alpha_1 T_t+\alpha_2 T_t^2+\alpha_3 T_t^3+u_t\\ =& \; 373580.083+14398.756 T_t+ 1537.276 T_t^2+167.452 T_t^3\\ =& \; 373580.083+14398.756 (97)+1537.276 (97)^2+167.452 (97)^3\\ =& \; 169'063,408 \end{aligned}\] La Cartera de Febrero de 2024 va ser de 169 mil millones.

Primer pronostico

Igual ahora con el modelo de tendencia cuadratica

  • Tomamos la ecuación del modelo (2) y se reemplaza en el último valor o periodo tendencial. P.e:

\[\begin{aligned} Cartera_t=& \; \alpha_0+\alpha_1 T_t+\alpha_2 T_t^2+u_t\\ =& \; 378280.9633+7663.445 T_t+ 3588.569 T_t^2\\ =& \; 378280.9633+7663.445 (97)+ 3588.569 (97)^2\\ =& \; 34'886.480 \end{aligned}\] La Cartera de Febrero de 2024 va ser de 34 mil millones.

Pronosticos con tendencia

datos <- read_excel("cartera.xlsx")
datos <- datos$Cartera
bsd2 <- ts(datos, start = c(2016, 1), frequency = 12)

# Tendencia lineal
m.lin <- tslm(bsd2 ~ trend)
p.lin <- forecast(m.lin, h=10)

# Tendencia Exponencial
m.exp <- tslm(bsd2 ~ trend, lambda = 0)
p.exp <- forecast(m.exp, h=10)

# Tendencia Cúbica
m.cub <- tslm(bsd2 ~ trend+trend^2+trend^3)
p.cub <- forecast(m.cub, h=10)
autoplot(bsd2) +
  autolayer(fitted(m.lin), series = "Lineal") +
  autolayer(fitted(m.exp), series="Exponencial") +
  autolayer(fitted(m.cub), series = "Cúbico") +
  autolayer(p.lin$mean, series = "Lineal") +
  autolayer(p.exp$mean, series="Exponencial") +
  autolayer(p.cub$mean, series="Cúbico") +
  xlab("Año") +  ylab("$ Miles de Millones") +
  ggtitle("Predicción simple") +
  guides(colour=guide_legend(title=" "))

Otros pronosticos

Media simple

Dado el comportamiento de la serie \(y_t\), se extrae la media o esperanza matemática y con eso se tiene el pronostico.

  • El pronostico de todos los valores futuros es igual a la media histórica de los datos \(y_1,\dots,y_T\).

  • La formula es: \(\hat{y}_{T+h|T} = \bar{y} = \frac{ (y_1+\dots+y_T)}{T}\)

Otros pronosticos

Ingenuo (Naive)

El pronostico es igual al último valor de la serie.

  • El pronostico de todos los valores futuros es igual al último dato \(y_T\).

  • La formula es: \(\hat{y}_{T+h|T} =y_T\)

Otros pronosticos

Método ingenuo estacional

Si se mira el comportamiento de la serie \(y_t\), El pronostico es similar al de la última season.

  • El pronostico de todos los valores futuros es igual al dato final de la estación del periodo \(y_{1s},\dots,y_{TS}\).

  • La formula es: \(\hat{y}_{T+h|T} =y_{T+h-m(k+1)}\), donde la parte de \(m=\) es el periodo estacional y \(k\) viene a ser un número producto de \(k=\frac{(h-1)}{m}\).

Otros pronosticos

Code 
# Grafico de pronósticos
autoplot(bsd2) +
  autolayer(meanf(bsd2, h=11), PI=FALSE, series="Promedio") +
  autolayer(naive(bsd2, h=11), PI=FALSE, series="Naïve") +
  autolayer(snaive(bsd2, h=11), PI=FALSE, series="Naïve Estacional") +
  ggtitle("Pronosticos de Cartera ") +
  xlab("Fecha") + ylab("$ en miles de millones") +
  guides(colour=guide_legend(title="Pronostico por"))

Box-cox 🤠

Modelo Box-cox

  • Es una manera de extraer una Tendencia no lineal

\[w_t = \left\{\begin{array}{ll} \log(y_t), & \quad \lambda = 0; \\ (y_t^\lambda-1)/\lambda , & \quad \lambda \ne 0. \end{array}\right.\]

  • Se sabe que si \(\lambda=1\): (No es necesario transformar)
  • \(\lambda=\frac12\): (Raíz cuadrada y adicional transformación lineal)
  • \(\lambda=0\): (Logaritmo neperiano)
  • \(\lambda=-1\): (Inversa)

Modelo Box-cox

Code 
autoplot(BoxCox(bsd2,lambda=1/3))

Modelo Box-cox

Adquiriendo el lambda correcto

  • Cuando queremos encontrar el lambda \((\lambda)\) optimizando el verdadero valor de acuerdo a la serie pero el sistema lo hace por nosotros 😃
Code 
(BoxCox.lambda(bsd2))
[1] -0.9999242
  • Valores pequeños del \(\lambda\) manifiestan intervalos mas grandes del pronostico

Pronostico con Box Cox

  • Si queremos obtener en valores reales el pronostivo, debemos reversar la transformación Box-Cox

\[ y_t = \left\{\begin{array}{ll} \exp(w_t), & \quad \lambda = 0; \\ (\lambda W_t+1)^{1/\lambda} , & \quad \lambda \ne 0. \end{array}\right.\]

Code 
niv1 <- snaive(bsd2, lambda=1/3)
autoplot(niv1)

Code 
niv2 <- snaive(bsd2, lambda=-0.999)
autoplot(niv2)

Gracias por su atención!!

Referencias

Cleveland, Robert B, William S Cleveland, Jean E McRae, and Irma Terpenning. 1990. “STL: A Seasonal-Trend Decomposition.” J. Off. Stat 6 (1): 3–73.
Hyndman, Rob J, and George Athanasopoulos. 2018. Forecasting: Principles and Practice. OTexts.

Universidad del Norte