PRIMERA PARTE

Las librerías y paquetes que utilizaremos para ejecutar todo nuestro análisis, son las siguientes:

library(ggplot2)
library(forecast)
library(tseries)
library(patchwork)
library(caret)
library(tidyverse)
library(robustbase)
library(cvTools)
require(season)
library(DAAG)

Ejercicio #1

Inicie con una contextualizacion de la serie: definición, frecuencia, periodo los datos que escogieron, fuente de los datos (solo si esta disponible)

Solución:

La serie temporal tomada para el presente análisis, corresponde al estudio de la popularidad de la empresa Argos según las búsquedas hechas en Google que contengan a esta empresa, para lo cual se toma dentro de un periodo de tiempo de 3 años comprendido entre 2017 y 2019 con una frecuencia de 5 datos por mes, dicha serie fue obtenida y descargada de la herramienta gratuita que ofrece google conocida como Google Trends, donde se ejecuta así:

Base_argos <-read.csv("C:/Users/Fabian Navarro/Google Drive/Taller_3_MMF/multiTimeline (1).csv")

Después de haber cargado esta base de datos, se puede obtener el siguiente gráfico de serie de tiempo:

popularidad <-ts(Base_argos$Argos, frequency = 5)
autoplot(popularidad, xlab = "Meses", main= "Popularidad de búsqueda")

Un análisis rápido y descriptivo de esta serie de tiempo seria lo siguiente:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   33.00   41.00   42.96   51.00  100.00

Ademas, el modelo se puede describir por una función ARIMA, siendo así :

## Series: popularidad 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     mean
##       0.8093  -0.4715  42.7106
## s.e.  0.0745   0.1062   2.9525
## 
## sigma^2 estimated as 189.2:  log likelihood=-632.98
## AIC=1273.95   AICc=1274.22   BIC=1286.18

En conclusión, con esto anterior se puede decir que nuestra serie de datos esta representada por un ARMA(1,1), es decir, donde su componente auto-regresiva es igual 1 y su componente de media móvil es también 1.

Ejercicio #2

Use la función \(stl()\) para determinar si existe componentes de tendencia y/o de estacionalidad.Reporte la gráfica y una conclusión preliminar sobre si existe tendencia global, local, y si existe componente estacional, de que periodo. Note que la tendencia global puede también ser una función periódica. En este caso habrían dos componentes estacionales, y la de mayor periodo, por ejemplo, anual, se puede combinar con una tendencia por ejemplo lineal

Solución:

En relación a la serie de tiempo anterior, si aplicamos la función \(stl()\) y graficamos su descomposición, se tendría lo siguiente:

Ahora bien, habiendo obtenido este gráfico, se puede presenciar que si existe estacionalidad y tendencia, tal que de manera preliminar la serie de datos presenta una tendencia global en la fecha aproximada del mes 24, es decir, en el mes de Diciembre del 2018, y de igual manera, su estacionalidad presenta patrones de recurrencia con un periodo aproximado de 5 meses.

Ejercicio #3

Decida cual modelo utilizar para la tendencia y cual para la estacionalidad. Escoja al menos dos modelos de la forma \(Y_{t}=T_{t}+S_{t}+ \epsilon_{t}\), siendo \(T_{t}\): lineal, cuadrática, cubica, exponencial-lineal, loess, filtro lineal; y \(S_{t}\): suma de variables indicadoras, suma de combinaciones seno y coseno. Es posible combinar estas componentes. Para definir la componente estacional hay que determinar primero el posible periodo o los posibles periodos. Al final, sugiero algunas consideraciones para la selección del modelo. Recuerde que no hay una regla que sirva para todas las series de manera que esta parte es también cuestión intuitiva. En cualquiera de los casos anteriores implemente la estrategia de validación cruzada escogiendo una parte de los datos para ajustar y otra para comparar los pronósticos.

Solución:

En este punto, nos enfocaremos en modelar la serie de tiempo anteriormente vista como una función de tendencia, aplicando así los siguientes modelos:

-> Modelo Lineal

mod_lin <-lm(popularidad~time(popularidad))
summary(mod_lin)
## 
## Call:
## lm(formula = popularidad ~ time(popularidad))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.834 -10.446  -2.894   7.648  58.169 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        45.6429     2.6266  17.377   <2e-16 ***
## time(popularidad)  -0.1615     0.1389  -1.163    0.247    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.77 on 155 degrees of freedom
## Multiple R-squared:  0.008651,   Adjusted R-squared:  0.002256 
## F-statistic: 1.353 on 1 and 155 DF,  p-value: 0.2466

En este primer modelo de tendencia, podemos ver que se presenta un error estándar de \(15.77\) y un error regresivo de \(2.6266\), con un intercepto en \(45.6429\), y ademas, una variabilidad igual a \(-0.1615\) lo que significa que con el pasar de los meses se esta perdiendo al rededor de \(0.1661\) búsquedas en Internet sobre la empresa Argos, lo cual hace que en teminos de tiempos largos la popularidad de búsqueda se vaya acabando.

-> Modelo Lineal Cuadrático

mod_cuad<-lm(popularidad~time(popularidad)+ I(time(popularidad)^2))
summary(mod_cuad)
## 
## Call:
## lm(formula = popularidad ~ time(popularidad) + I(time(popularidad)^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.317 -10.425  -3.126   7.480  57.867 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            43.880688   4.237729  10.355   <2e-16 ***
## time(popularidad)       0.141000   0.586803   0.240    0.810    
## I(time(popularidad)^2) -0.009112   0.017170  -0.531    0.596    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.81 on 154 degrees of freedom
## Multiple R-squared:  0.01046,    Adjusted R-squared:  -0.00239 
## F-statistic: 0.814 on 2 and 154 DF,  p-value: 0.445

Para este otro modelo de tendencia, podemos ver que se presenta un error estándar de \(15.81\) y un error regresivo de \(4.2377\), con un intercepto aproximado de \(43.8807\), y que ademas, tiene una variabilidad en su componente cuadrática igual a \(-0.0091\) lo que significa que con el pasar de los meses se esta perdiendo al rededor de \(0.0091\) búsquedas en Internet sobre la empresa Argos, lo cual hace que en teminos de tiempos largos la popularidad de búsqueda se vaya acabando.

-> Modelo Lineal Cubico

mod_cub<-lm(popularidad~time(popularidad)+I(time(popularidad)^2)+I(time(popularidad)^3))
summary(mod_cub)
## 
## Call:
## lm(formula = popularidad ~ time(popularidad) + I(time(popularidad)^2) + 
##     I(time(popularidad)^3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.941  -9.938  -1.292   7.698  52.458 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            60.444136   5.981759  10.105  < 2e-16 ***
## time(popularidad)      -5.163709   1.513653  -3.411 0.000827 ***
## I(time(popularidad)^2)  0.380067   0.104388   3.641 0.000371 ***
## I(time(popularidad)^3) -0.007815   0.002070  -3.776 0.000228 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.17 on 153 degrees of freedom
## Multiple R-squared:  0.0948, Adjusted R-squared:  0.07705 
## F-statistic: 5.341 on 3 and 153 DF,  p-value: 0.001587

Para este otro modelo de tendencia, podemos ver que se presenta un error estándar de \(15.17\) y un error regresivo de \(5.9817\), con un intercepto aproximado de \(60.4441\), y que adicionalmente, tiene una variabilidad en su componente cubica igual a \(-0.0078\) lo que significa que con el pasar de los meses se esta perdiendo al rededor de \(0.0078\) búsquedas en Internet sobre la empresa Argos, lo cual hace que en teminos de tiempos largos la popularidad de búsqueda se vaya acabando.

-> Modelo de regresión Loess

mod_loess <- loess(popularidad~time(popularidad))
summary(mod_loess)
## Call:
## loess(formula = popularidad ~ time(popularidad))
## 
## Number of Observations: 157 
## Equivalent Number of Parameters: 4.37 
## Residual Standard Error: 15.06 
## Trace of smoother matrix: 4.76  (exact)
## 
## Control settings:
##   span     :  0.75 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE

Como ultimo, se tiene que este modelo de tendencia, presenta un error estándar de \(15.06\) y una traza de suavizamiento del \(4.76\).

(Importante) Habiendo ya analizado estos modelos de tendencia, según los resultados obtenidos, consideramos que los dos mejores modelos que se acercan a nuestra serie de tiempo son los modelos de: “Regresión Loess” y “Lineal Cubico”, ya que estos presentan menor error estándar y se ajustan adecuadamente para así realizar un suavizamiento optimo.

Por otro lado, si nos enfocamos en modelar esta serie de tiempo como una función de estacionalidad, se puede aplicar el siguiente modelo:

-> Modelo para estacionalidad

It2 = fourier(popularidad,2)
mod_st = lm(popularidad ~ time(popularidad) + It2)
summary(mod_st)
## 
## Call:
## lm(formula = popularidad ~ time(popularidad) + It2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.011 -10.367  -2.903   6.903  58.160 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       45.63210    2.65794  17.168   <2e-16 ***
## time(popularidad) -0.16114    0.14053  -1.147    0.253    
## It2S1-5           -0.01438    1.79859  -0.008    0.994    
## It2C1-5           -1.06311    1.80434  -0.589    0.557    
## It2S2-5           -0.38832    1.79869  -0.216    0.829    
## It2C2-5           -0.11999    1.80434  -0.067    0.947    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.96 on 151 degrees of freedom
## Multiple R-squared:  0.01127,    Adjusted R-squared:  -0.02146 
## F-statistic: 0.3444 on 5 and 151 DF,  p-value: 0.8852

Ejercicio #4

Estime al menos dos modelos. Si incluyo el modelo exponencial debe estimar también el \(log-lineal\) como auxiliar. Reporte las tablas de parámetros estimados, estadísticos \(t\) y \(F\), valores \(p\). Finalmente, concluya si los modelos ajustan

Solución:

En este siguiente ejercicio, estimaremos otros dos modeles no vistos en el anterior punto, pero que incluso tiene relación en algunos cálculos de este previo punto, por tal razón, se tiene los siguientes modelos:

-> Modelo \(log-lineal\)

mod_log_lin<-lm(log(popularidad)~time(popularidad))
summary(mod_log_lin)
## 
## Call:
## lm(formula = log(popularidad) ~ time(popularidad))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07849 -0.21759 -0.00624  0.21173  0.95419 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.799226   0.060739  62.550   <2e-16 ***
## time(popularidad) -0.006282   0.003211  -1.956   0.0523 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3647 on 155 degrees of freedom
## Multiple R-squared:  0.02409,    Adjusted R-squared:  0.01779 
## F-statistic: 3.826 on 1 and 155 DF,  p-value: 0.05226

En este modelo se puede verificar que su valor \(R^{2}\) normal es igual a \(0.0241\) para un ajuste optimo del modelo, de igual manera, también se visualiza que su \(Valor-p\) es igual a \(0.0523\) implicando que sus datos son no representativo para el modelo, por ende, no se rechaza la hipótesis nula y por ultimo, su valor \(F\) (Fisher) es igual a \(3.826\)

-> Modelo \(Exp-lineal\)

Ds= data.frame(popularidad,time(popularidad))
beta0= mod_log_lin$coefficient[1]
beta1= mod_log_lin$coefficient[2]
mod_exp= nls(popularidad~exp(beta0+beta1*time(popularidad)),data=Ds,start=list(beta0=beta0,beta1=beta1))

exp_ln <- lm(mod_exp)
summary(exp_ln)
## 
## Call:
## lm(formula = mod_exp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.797 -10.404  -2.921   7.604  58.172 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             17.5144    22.1532   0.791    0.430
## exp(beta0 + beta1 * time(popularidad))   0.6313     0.5487   1.151    0.252
## 
## Residual standard error: 15.77 on 155 degrees of freedom
## Multiple R-squared:  0.008468,   Adjusted R-squared:  0.002071 
## F-statistic: 1.324 on 1 and 155 DF,  p-value: 0.2517

En este ultimo modelo se puede verificar que su valor \(R^{2}\) normal es igual a \(0.0085\) para un ajuste optimo del modelo, de igual forma, también se puede visualizar que su \(Valor-p\) es igual a \(0.2517\) que implicando que sus datos son no representativo para el modelo, por ende, no se determina rechazar la hipótesis nula y por ultimo, su valor \(F\) (Fisher) es igual a \(1.324\)

(Importante) A partir de estas cifras analizadas de los dos modelos anteriores, se puede decir que el modelo \(log-lineal\) tiene el valor \(R^{2}\) normal mas “próximo” a \(1\), implicando que este modelo se ajuste mejor a los datos de la serie de tiempo.

Ejercicio #5

Reporte AIC, BIC y R cuadrado ajustado para los modelos. Escoja uno de los modelos con base en estos criterios

Solución:

Nuestro reporte de AIC, BIC y R, para los modelos lineal, lineal cuadrático, lineal cubico y \(exp-lineal\), se reporta tal como se ve la siguiente tabla:

AIC.tot= c(AIC(mod_lin), AIC(mod_cuad),AIC(mod_cub), AIC(exp_ln))
BIC.tot= c(AIC(mod_lin,k=log(5)),AIC(mod_cuad,k=log(5)),AIC(mod_cub,k=log(5)),AIC(exp_ln,k=log(5)))

R.tot=c(summary(mod_lin)$r.squared,summary(mod_cuad)$r.squared,summary(mod_cub)$r.squared,summary(exp_ln)$r.squared)

medidas1 = rbind(AIC.tot, BIC.tot, R.tot)
colnames(medidas1) = c("lin","cuad","cub","exp_ln")
rownames(medidas1) = c("AIC","BIC","R-squared")
medidas1
##                    lin         cuad          cub       exp_ln
## AIC       1.315629e+03 1.317342e+03 1.305356e+03 1.315658e+03
## BIC       1.314457e+03 1.315779e+03 1.303403e+03 1.314486e+03
## R-squared 8.651330e-03 1.046089e-02 9.479632e-02 8.468244e-03
#knitr::kable((medidas1))

Ejercicio #6

Calcule los pronósticos para la validación cruzada con los modelos escogidos. Calcule el MAPE para cada modelo. Reporte la gráfica de los pronósticos y de los datos observados, para el periodo de comparación. ¿Cual modelo es el escogido para pronosticar?

Solución:

En este punto realizaremos un pronostico por medio de la validación cruzada “k-fold” haciendo uso del paquete forecast, obteniendo así lo siguiente:

far2 <- function(popularidad, h=5){forecast(Arima(popularidad, order=c(1,0,1)), h=5)}
modelcv <- CVar(popularidad,k=5, lambda=0.152, fit= far2)
print(modelcv)
## Series: popularidad 
## Call:   CVar(y = popularidad, k = 5, lambda = 0.152, fit = far2)
## 
## 5-fold cross-validation
##                 Mean         SD
## ME        -1.4070794  3.6217614
## RMSE      15.9483643  3.3305938
## MAE       12.5459562  2.8041179
## MPE       -6.7333106 10.0998592
## MAPE      31.3107740  7.8427803
## ACF1       0.2164175  0.2042129
## Theil's U  1.2973416  0.2105491
## 
## p-value of Ljung-Box test of residuals is  0.741559 
## if this value is significant (<0.05),
## the result of the cross-validation should not be used
## as the model is underfitting the data.

Ejercicio #7

Aplique el método de Holt-Winters y calcule los valores suavizados y los pronósticos. Compare los resultados de este método con los de los dos modelos de regresión escogidos. Reporte las gráficas de los valores suavizados y de la serie superpuestos. Reporte el MAPE y las otras medidas de error de los pronósticos, calculadas con la función accuracy

Solución:

En este ejerció, debemos aplicar el método Holt-Winters, tal que esto nos indique como son los ajustes idóneos para una optima predicción, estos se visualizan de la siguiente manera:

Holt_W = HoltWinters(popularidad)
pre_Holt = predict(Holt_W, 8, prediction.interval = TRUE)
print(pre_Holt)
## Time Series:
## Start = c(32, 3) 
## End = c(33, 5) 
## Frequency = 5 
##           fit      upr        lwr
## 32.4 33.41991 62.88095  3.9588741
## 32.6 32.30279 63.76797  0.8376024
## 32.8 33.96908 67.31819  0.6199789
## 33.0 35.86885 71.00100  0.7366956
## 33.2 33.01990 69.84887 -3.8090742
## 33.4 35.33420 74.78472 -4.1163257
## 33.6 34.21707 75.18595 -6.7518041
## 33.8 35.88337 78.31630 -6.5495633

Ahora bien, habiendo ya realizado este análisis de valores para la predicción, si graficamos su aproxmiacion, podríamos ver esto:

plot(time(popularidad),popularidad,type='l', col='red',ylab='log(y)')
lines(c(time(popularidad),seq(1+T,8+T)),c(popularidad[1:5],Holt_W$fitted[,1],pre_Holt[,1]),col = 'blue')
legend("topleft", legend = c("Serie Original", "Holt-Winters"),lwd = 3, col = c("red", "blue"))

Por otro lado, si realizamos el análisis MAPE, tendríamos lo siguiente:

Acu=accuracy(popularidad[1:24],pre_Holt[,1])
summary(Acu)
##        ME           RMSE           MAE          MPE              MAPE      
##  Min.   :-16   Min.   :18.6   Min.   :16   Min.   :-46.59   Min.   :46.59  
##  1st Qu.:-16   1st Qu.:18.6   1st Qu.:16   1st Qu.:-46.59   1st Qu.:46.59  
##  Median :-16   Median :18.6   Median :16   Median :-46.59   Median :46.59  
##  Mean   :-16   Mean   :18.6   Mean   :16   Mean   :-46.59   Mean   :46.59  
##  3rd Qu.:-16   3rd Qu.:18.6   3rd Qu.:16   3rd Qu.:-46.59   3rd Qu.:46.59  
##  Max.   :-16   Max.   :18.6   Max.   :16   Max.   :-46.59   Max.   :46.59  
##       ACF1           Theil's U    
##  Min.   :-0.3229   Min.   :10.48  
##  1st Qu.:-0.3229   1st Qu.:10.48  
##  Median :-0.3229   Median :10.48  
##  Mean   :-0.3229   Mean   :10.48  
##  3rd Qu.:-0.3229   3rd Qu.:10.48  
##  Max.   :-0.3229   Max.   :10.48

Ejercicio #8

Calcule la media móvil doble tipo Henderson de longitud \(13\), aplicada a la serie, denotada por m.r y definida como sigue: \[\begin{align*} &w13 = c(−0,019,−0,028,0,0,0,066,0,147,0,214,0,240,0,214,0,147,0,066,0,0,−0,028,−0,019)\\ &m.r=na.omit(filter(y, w13,``conv",2, F,NULL)) \end{align*}\]

Reporte las gráficas de la tendencia \(T_{t}\) con el modelo de regresión escogido y de la media móvil \(m.r\), superpuestas. Comente sobre la diferencia. ¿Cual captura mejor la tendencia de la serie?

Defina el vector de incrementos porcentuales de \(m.r\) como \(icp.r = (m.r[2 : n] − m.r[1 : (n − 1)])/m.r[1 : (n − 1)]\), donde \(n = length(m.r)\). Calcule el promedio de estos valores icp.r. En caso de que la tendencia sea monótona (creciente o decreciente) ¿como puede interpretarse esta media?

Solución:

En este ultimo ejercicio, si ejecutamos una media móvil de longitud \(13\) al modelo de Holt-Winters, se tendría la siguiente gráficas:

#w13 <- c(-0.019, -0.028, 0.0, 0.066, 0.147, 0.214, 0.240, 0.214, 0.147, 0.066, 0.0, -0.028, -0.019)
#m.r <- na.omit(filter(popularidad, w13, "conv", 2, F,NULL))
#summary(m.r)

Habiendo realizado el código anterior para la media móvil de longitud \(13\) dentro del modelo de Holt-Winters, se tiene como análisis preliminar, que este modelo presenta una media equivalente a \(42.91\), y aun así si se verifica como es el comportamiento de este método gráficamente, se tiene lo siguiente :

#plot(m.r,col='red',lwd=2)
#lines(mod_loess)
#legend("topleft", legend = c("m.r", "loess"), lwd = 3, col = c("red", "black"))

Ahora bien, respecto a las gráfica anterior se puede percibir que el modelo de Holt-Winters (m.r), presenta una mejor aproximación o ajuste a la tendencia de la serie de tiempo original. Consiguientemente, si definimos el vector de incrementos porcentuales de \(m.r\), se tendría:

#n = length(m.r)
#icp.r = (m.r[2 : n] - m.r[1 : (n - 1)])/m.r[1 : (n - 1)]
#plot(icp.r)
#summary(icp.r)

Por ultimo, en referencia al resumen anterior y la gráfica de tendencia de la serie original, se puede decir que en promedio cada ciclo de nuestra serie esta incrementando en un \(0.0003948\).

Conclusiones generales del taller

En este trabajo, se evidencio que los modelos no se ajustaban optimamente a la serie original, debido a que los valores \(R^{2}\) no se acercaban lo suficientemente a \(1\), pero, como análisis exploratorio del taller se llego a que la serie suavizada logaritmicamente presenta una disminución de los errores, brindado valores totalmente distintos.