Obtener Datos

El fichero Estadísticas.xls contiene la serie temporal de jubilados en España en los años 2005 al presente.

Se pide:

  1. Leer el fichero

  2. Convertir en objeto ts la serie temporal de jubilados. Representar y analizar la serie.

setwd("/cloud/project/Modulo 6/Examenes")
datos=read.csv("Estadisticas.csv",sep=",")
tail(datos)
##              Fecha Jubilados
## 227 Noviembre 2023   6297621
## 228 Diciembre 2023   6307217
## 229     Enero 2024   6328280
## 230   Febrero 2024   6334392
## 231     Marzo 2024   6341876
## 232     Abril 2024   6348237
jubilados=ts(datos[,2],frequency = 12,star=c(2005,1))
summary(jubilados)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 4582093 4939833 5430886 5416767 5899444 6348237
plot(jubilados)

La serie temporal de jubilados presenta una tendencia creciente a lo largo del tiempo, no es estacionaria en media, y no se observa estacionalidad anual.

  1. Obtener una primera diferencia regular.Representar y analizar la serie.
jubilados_d=diff(jubilados,1)
summary(jubilados_d)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -20179    5160    8314    7607   10870   21063
plot(jubilados_d)

La primera diferencia hace estacionaria la serie en media y varianza, hay probablemente un outlier en mayo de 2020, y se observa una estructura cíclica.

Modelo ARIMA

Estime un modelo ARIMA para la serie de jubilados;

  1. Estudie la estacionariedad de la serie
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
adf.test(jubilados, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  jubilados
## Dickey-Fuller = -2.8485, Lag order = 6, p-value = 0.2193
## alternative hypothesis: stationary
ndiffs(jubilados)
## [1] 1
nsdiffs(jubilados)
## [1] 1

El test de Dickey-Fuller no rechaza la no estacionariedad de la serie original, de hecho las pruebas sobre estacionariedad aconsejan diferenciarla.

  1. Realice las transformaciones necesarias para hacerla estacionaria
adf.test(jubilados_d, alternative = "stationary")
## Warning in adf.test(jubilados_d, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  jubilados_d
## Dickey-Fuller = -5.8991, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
ndiffs(jubilados_d)
## [1] 0
nsdiffs(jubilados_d)
## [1] 0

La estacionariedad de la serie, se estudia con el test de Dickey-Fuller aumentado, cuyo p-valor rechaza la hipótesis nula de no estacionariedad de la serie. La serie con una diferencia regular es una serie estacionaria, y que por tanto no precisa de ninguna diferenciación adicional.

  1. Identifique un modelo ARIMA
Acf(jubilados_d)

Pacf(jubilados_d)

En la parte regular podemos comprobar que las funciones de correlación simple el primer retardo y el segundo son significativos y en la función de autocorrelación parcial podemos considerar solo el primer retardo significativos. Esto indica que quizás tengamos que tener en cuenta modelos con MA1 y AR1 y/o AR2

Por otro lado, en la parte estacional vemos tanto en la acf como en pacf el retardo 12 es significativo y el 24 en la paf.

Posibles modelos: \((2,1,1) (1,0,1)_{12}\), \((1,1,1) (1,0,1)_{12}\), \((1,1,1) (1,0,2)_{12}\)

  1. Estime el modelo ARIMA
mod1=arima(jubilados, order = c(2,1,1), seasonal = list(order = c(1,0,1)))
mod1 
## 
## Call:
## arima(x = jubilados, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 1)))
## 
## Coefficients:
##          ar1     ar2      ma1    sar1     sma1
##       0.8231  0.0792  -0.6752  0.9426  -0.6150
## s.e.  0.1478  0.0970   0.1350  0.0337   0.0985
## 
## sigma^2 estimated as 16610051:  log likelihood = -2253.99,  aic = 4519.98
plot(mod1)

T<-mod1$coef/sqrt(diag(mod1$var.coef))
T
##        ar1        ar2        ma1       sar1       sma1 
##  5.5685213  0.8157821 -5.0002665 27.9763842 -6.2448480

Se obtienen coeficiente no significativos en ar2

Raices dentro del circulo unitario

mod2=arima(jubilados, order = c(1,1,1), seasonal = list(order = c(1,0,1)))
mod2 
## 
## Call:
## arima(x = jubilados, order = c(1, 1, 1), seasonal = list(order = c(1, 0, 1)))
## 
## Coefficients:
##          ar1      ma1    sar1     sma1
##       0.9302  -0.7472  0.9398  -0.6161
## s.e.  0.0435   0.0781  0.0357   0.1028
## 
## sigma^2 estimated as 16685832:  log likelihood = -2254.33,  aic = 4518.66
plot(mod2)

Se obtienen coeficientes significativos Raices dentro del circulo unitario

T<-mod2$coef/sqrt(diag(mod2$var.coef))
T
##       ar1       ma1      sar1      sma1 
## 21.372741 -9.570602 26.288869 -5.993015

coeficientes significativos

mod3=arima(jubilados, order = c(1,1,1), seasonal = list(order = c(1,0,2)))
mod3 
## 
## Call:
## arima(x = jubilados, order = c(1, 1, 1), seasonal = list(order = c(1, 0, 2)))
## 
## Coefficients:
##          ar1      ma1    sar1     sma1     sma2
##       0.9181  -0.7260  0.9670  -0.6035  -0.1181
## s.e.  0.0505   0.0854  0.0241   0.0789   0.0756
## 
## sigma^2 estimated as 16423803:  log likelihood = -2253.23,  aic = 4518.47
plot(mod3)

Raices dentro del circulo unitario

T<-mod3$coef/sqrt(diag(mod3$var.coef))
T
##       ar1       ma1      sar1      sma1      sma2 
## 18.171171 -8.496346 40.119510 -7.647324 -1.563399

Se obtiene un coeficiente no significativo en sma2

Buscamos el modelo que propone autoarima

mod4=auto.arima(jubilados)
mod4
## Series: jubilados 
## ARIMA(1,1,2)(0,1,2)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sma1     sma2
##       0.7901  -0.6555  0.0754  -0.6081  -0.1382
## s.e.  0.1213   0.1419  0.0707   0.0706   0.0703
## 
## sigma^2 = 18122073:  log likelihood = -2134.81
## AIC=4281.63   AICc=4282.02   BIC=4301.96
T<-mod4$coef/sqrt(diag(mod4$var.coef))
T
##       ar1       ma1       ma2      sma1      sma2 
##  6.514175 -4.618930  1.067455 -8.615456 -1.967203

Se obtienen coeficientes no significativos ma2 y sma2

Me quedo con el modelo: \((1,1,1) (1,0,1)_{12}\)

  1. Evalue los resultados

Analisis de los residuos

tsdisplay(residuals(mod2), lag.max=10, main="Residuos del modelo")

Acf(residuals(mod2))

Pacf(residuals(mod2))

Los residuos no tienen estructura.

Test de normalidad

hist(mod2$residuals)

jarque.bera.test(mod3$residuals)
## 
##  Jarque Bera Test
## 
## data:  mod3$residuals
## X-squared = 66.656, df = 2, p-value = 3.331e-15
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.96717, p-value = 3.398e-05

Los residuos no son normales, el histograma muestra asimetria.

  1. Realice una estimación para lo que queda de año
pronostico <- forecast(mod2,h=8)
plot(pronostico)

print(pronostico)
##          Point Forecast   Lo 80   Hi 80   Lo 95   Hi 95
## May 2024        6348732 6343497 6353967 6340726 6356739
## Jun 2024        6361135 6353025 6369244 6348733 6373537
## Jul 2024        6371596 6360828 6382364 6355128 6388064
## Aug 2024        6380254 6366891 6393618 6359817 6400692
## Sep 2024        6387260 6371322 6403197 6362885 6411634
## Oct 2024        6398616 6380111 6417120 6370315 6426916
## Nov 2024        6412242 6391175 6433309 6380023 6444462
## Dec 2024        6424052 6400427 6447677 6387921 6460183

Desconposicion temporal

  1. Realice la descomposición temporal de la serie de jubilados con los metodos, stl, RJDemetra y Prophet.
library(forecast)
mod5=stl(jubilados,"per")
plot(stl(jubilados,"per"))

str(mod5$time.series)
##  Time-Series [1:232, 1:3] from 2005 to 2024: 8150.2 2069.9 14.1 -1480.1 -8557.6 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:3] "seasonal" "trend" "remainder"
head(mod5$time.series)
##             seasonal   trend remainder
## Jan 2005  8150.24513 4579845  3120.632
## Feb 2005  2069.90097 4586019 -5996.383
## Mar 2005    14.05573 4592194 -6512.897
## Apr 2005 -1480.09984 4598273 -4892.631
## May 2005 -8557.59578 4604352  7158.976
## Jun 2005 -5334.05266 4610396  5401.656
fcs_stl = forecast(mstl(jubilados), h = 8)
autoplot(fcs_stl)

RJ Demetra

library(RJDemetra)
spec_ts <- tramoseats_spec(spec = "RSA3")
mod6 <- tramoseats(jubilados, spec = spec_ts)
mod6
## RegARIMA
## y = regression model + arima (3, 1, 1, 1, 0, 1)
## Log-transformation: yes
## Coefficients:
##           Estimate Std. Error
## Phi(1)    -0.46833      0.218
## Phi(2)    -0.15270      0.073
## Phi(3)    -0.09785      0.092
## BPhi(1)   -0.84885      0.052
## Theta(1)  -0.48412      0.216
## BTheta(1) -0.42930      0.084
## 
##              Estimate Std. Error
## Mean         0.001425      0.000
## LS (5-2020) -0.003610      0.001
## LS (2-2005) -0.002657      0.001
## 
## 
## Residual standard error: 0.000683 on 221 degrees of freedom
## Log likelihood =  1352, aic =  4477 aicc =  4478, bic(corrected for length) = -14.37
## 
## 
## 
## Decomposition
## Model
## AR :  1 - 0.468330 B - 0.152697 B^2 - 0.097848 B^3 - 0.848847 B^12 + 0.397540 B^13 + 0.129616 B^14 + 0.083058 B^15 
## D :  1 - B 
## MA :  1 - 0.484119 B - 0.429300 B^12 + 0.207832 B^13 
## 
## 
## SA
## AR :  1 - 1.793929 B + 0.796540 B^2 
## D :  1 - B 
## MA :  1 - 1.738404 B + 0.907207 B^2 - 0.145047 B^3 
## Innovation variance:  0.5799575 
## 
## Trend
## AR :  1 - 1.793929 B + 0.796540 B^2 
## D :  1 - B 
## MA :  1 - 0.481394 B - 0.969353 B^2 + 0.512040 B^3 
## Innovation variance:  0.08711419 
## 
## Seasonal
## AR :  1 + 1.325600 B + 1.428795 B^2 + 1.409415 B^3 + 1.390299 B^4 + 1.371441 B^5 + 1.352840 B^6 + 1.334490 B^7 + 1.316390 B^8 + 1.298535 B^9 + 1.280922 B^10 + 1.263548 B^11 + 0.397563 B^12 + 0.104273 B^13 
## MA :  1 + 1.666981 B + 2.063316 B^2 + 2.485695 B^3 + 2.716978 B^4 + 2.749006 B^5 + 2.579515 B^6 + 2.199943 B^7 + 1.658775 B^8 + 0.947749 B^9 + 0.339925 B^10 - 0.074262 B^11 - 0.141105 B^12 - 0.092562 B^13 
## Innovation variance:  0.1189709 
## 
## Irregular
## Innovation variance:  0.161608 
## 
## 
## 
## Final
## Last observed values
##                y      sa       t         s         i
## May 2023 6223630 6230140 6231359 0.9989550 0.9998045
## Jun 2023 6237879 6241343 6241690 0.9994449 0.9999445
## Jul 2023 6250037 6252666 6253128 0.9995796 0.9999261
## Aug 2023 6262124 6265635 6264846 0.9994397 1.0001260
## Sep 2023 6269845 6276043 6276067 0.9990124 0.9999963
## Oct 2023 6283169 6287057 6286661 0.9993815 1.0000631
## Nov 2023 6297621 6297110 6296578 1.0000812 1.0000845
## Dec 2023 6307217 6304623 6306585 1.0004114 0.9996890
## Jan 2024 6328280 6318636 6317314 1.0015263 1.0002093
## Feb 2024 6334392 6327747 6327617 1.0010501 1.0000205
## Mar 2024 6341876 6337141 6337167 1.0007471 0.9999959
## Apr 2024 6348237 6346434 6346655 1.0002842 0.9999651
## 
## Forecasts:
##              y_f    sa_f     t_f       s_f i_f
## May 2024 6350361 6356277 6356277 0.9990693   1
## Jun 2024 6362623 6365973 6365973 0.9994737   1
## Jul 2024 6373057 6375669 6375669 0.9995903   1
## Aug 2024 6382294 6385367 6385367 0.9995188   1
## Sep 2024 6389790 6395067 6395067 0.9991748   1
## Oct 2024 6401397 6404771 6404771 0.9994733   1
## Nov 2024 6414877 6414478 6414478 1.0000621   1
## Dec 2024 6426301 6424191 6424191 1.0003285   1
## Jan 2025 6442331 6433908 6433908 1.0013091   1
## Feb 2025 6449439 6443632 6443632 1.0009013   1
## Mar 2025 6457597 6453361 6453361 1.0006564   1
## Apr 2025 6464624 6463097 6463097 1.0002363   1
## 
## 
## Diagnostics
## Relative contribution of the components to the stationary
## portion of the variance in the original series,
## after the removal of the long term trend
##  Trend computed by Hodrick-Prescott filter (cycle length = 8.0 years)
##            Component
##  Cycle        52.967
##  Seasonal      6.120
##  Irregular     0.301
##  TD & Hol.     0.000
##  Others       22.531
##  Total        81.919
## 
## Combined test in the entire series
##  Non parametric tests for stable seasonality
##                                                           P.value
##    Kruskall-Wallis test                                          0
##    Test for the presence of seasonality assuming stability       0
##    Evolutive seasonality test                                    0
##  
##  Identifiable seasonality present
## 
## Residual seasonality tests
##                                       P.value
##  qs test on sa                          1.000
##  qs test on i                           1.000
##  f-test on sa (seasonal dummies)        1.000
##  f-test on i (seasonal dummies)         1.000
##  Residual seasonality (entire series)   1.000
##  Residual seasonality (last 3 years)    0.999
##  f-test on sa (td)                      0.471
##  f-test on i (td)                       0.112
## 
## 
## Additional output variables

RJ Demetra valora como outlier los datos del 5 del 2020 y 2 del 2005. Realiza transformacion logaritmica, presenta normalidad en los residuos:

jarque.bera.test(mod6$regarima$residuals)
## 
##  Jarque Bera Test
## 
## data:  mod6$regarima$residuals
## X-squared = 4.2156, df = 2, p-value = 0.1215
shapiro.test(mod6$regarima$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod6$regarima$residuals
## W = 0.98379, p-value = 0.0097

Prophet

library(prophet)
## Loading required package: Rcpp
## Loading required package: rlang
ds = seq(as.Date('2005-01-01'), as.Date('2024-04-01'), by = 'm')
jubilados1=data.frame(ds,jubilados)
colnames(jubilados1)=c("ds","y")
head(jubilados1)
##           ds       y
## 1 2005-01-01 4591116
## 2 2005-02-01 4582093
## 3 2005-03-01 4585695
## 4 2005-04-01 4591900
## 5 2005-05-01 4602953
## 6 2005-06-01 4610464
m<- prophet(jubilados1)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
future <- make_future_dataframe(m, periods = 8,freq = "months", include_history = TRUE)
head(future)
##           ds
## 1 2005-01-01
## 2 2005-02-01
## 3 2005-03-01
## 4 2005-04-01
## 5 2005-05-01
## 6 2005-06-01
fcs_php <- predict(m, future)
head(fcs_php)
##           ds   trend additive_terms additive_terms_lower additive_terms_upper
## 1 2005-01-01 4583295      2530.1280            2530.1280            2530.1280
## 2 2005-02-01 4589556     -7639.3695           -7639.3695           -7639.3695
## 3 2005-03-01 4595211       176.0727             176.0727             176.0727
## 4 2005-04-01 4601472     -1266.0237           -1266.0237           -1266.0237
## 5 2005-05-01 4607531     -9058.6774           -9058.6774           -9058.6774
## 6 2005-06-01 4613792     -7654.8439           -7654.8439           -7654.8439
##       yearly yearly_lower yearly_upper multiplicative_terms
## 1  2530.1280    2530.1280    2530.1280                    0
## 2 -7639.3695   -7639.3695   -7639.3695                    0
## 3   176.0727     176.0727     176.0727                    0
## 4 -1266.0237   -1266.0237   -1266.0237                    0
## 5 -9058.6774   -9058.6774   -9058.6774                    0
## 6 -7654.8439   -7654.8439   -7654.8439                    0
##   multiplicative_terms_lower multiplicative_terms_upper yhat_lower yhat_upper
## 1                          0                          0    4576496    4594460
## 2                          0                          0    4573381    4590621
## 3                          0                          0    4587361    4604526
## 4                          0                          0    4591473    4608955
## 5                          0                          0    4589738    4606882
## 6                          0                          0    4597908    4615553
##   trend_lower trend_upper    yhat
## 1     4583295     4583295 4585826
## 2     4589556     4589556 4581917
## 3     4595211     4595211 4595387
## 4     4601472     4601472 4600206
## 5     4607531     4607531 4598472
## 6     4613792     4613792 4606137
plot(m, fcs_php)

prophet_plot_components(m, fcs_php)

  1. Presente en una tabla los pronosticos realizados con los procedimientos que ha seguido, elija uno de ellos y justifique su elección.
time = tail(future$ds,8)
pronostico_stl_df <- data.frame(stl = fcs_stl$mean)
pronostico_php_df <- data.frame(prophet =tail( fcs_php$yhat,8))
pronostico_rj_df <- data.frame(rj = head(mod6$final$forecasts[,1],8))
pronosticos <- data.frame(time,pronostico_stl_df$stl,pronostico_php_df$prophet,pronostico_rj_df)
colnames(pronosticos) <- c("time", "stl", "prophet", "rj")
pronosticos
##         time     stl prophet      rj
## 1 2024-05-01 6344263 6327075 6350361
## 2 2024-06-01 6358265 6335043 6362623
## 3 2024-07-01 6370519 6345194 6373057
## 4 2024-08-01 6381015 6356235 6382294
## 5 2024-09-01 6389787 6361816 6389790
## 6 2024-10-01 6403673 6374401 6401397
## 7 2024-11-01 6420259 6389140 6414877
## 8 2024-12-01 6435579 6400368 6426301

Obtengo errores absolutos entre la serie original y los tres pornosticos:

prophet_error=sum(abs(jubilados-c(fcs_php$yhat)[1:232]))
prophet_error
## [1] 1130396
stl_error=sum(abs(jubilados-fcs_stl$fitted))
stl_error
## [1] 579160
rj_error=sum(abs(jubilados-mod6$final$series[,3]))
rj_error
## [1] 840748.6

La serie que presenta un error absoluto mas bajo es stl. Elegimos por tanto stl.

Suavizados

Vamos a suavizar la serie de diferencias mensuales de jubilados utilizando:

  1. el método de medias móviles. Justifique en tamaño de la media móvil elegido.
#Medias Móviles
#Alisado MA centrados
alisado1 = filter(jubilados_d, rep(1,12)/12, side=2)
#Alisado MA no centrados
alisado2 = filter(jubilados_d, rep(1,3)/3, side=1)
#Hacemos la gráfica para comparar
plot.ts(jubilados_d,main="Jubilados mes",xlab="Tiempo",ylab="personas")
lines(alisado1, col="red")
lines(alisado2, col="blue")
legend("bottomleft", c("Original", "Media móvil centrada 12",
   "Media móvil centrada 6"),
   lwd=c(1,2,2), col=c("black", "red", "blue"))
grid()

La media movil centrada de 12 terminos es menos suavizada que la se 6 terminos. Se elige la de 6 términos.

  1. los polinomios locales (kernel). Justifique el ancho de ventanda elegido.
library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
n=length(jubilados_d)
plot(1:n,jubilados_d, ylab="Velocidad", xlab="Distancia", main="Regresión kernel con diferentes parámetros de suavizado")
lines(ksmooth(1:n,jubilados_d, "normal", bandwidth = 2), col = 2)
lines(ksmooth(1:n,jubilados_d, "normal", bandwidth = 5), col = 3)
lines(ksmooth(1:n,jubilados_d, "normal", bandwidth = dpill(1:n,jubilados_d)), col = 4)

Se representan diferentes anchos de ventana, inclido el que minimiza el error cuadrático medio (dpill),

  1. los splines. Justifique el ancho de ventana elegido. Pronostique las personas que se jubilaran en lo que queda del año
plot(1:n,jubilados_d)
lines(smooth.spline(1:n,jubilados_d,spar=0.1),col="blue")
lines(smooth.spline(1:n,jubilados_d),col="red")

spl=smooth.spline(1:n,jubilados_d)
spl
## Call:
## smooth.spline(x = 1:n, y = jubilados_d)
## 
## Smoothing Parameter  spar= 0.3781272  lambda= 8.316674e-07 (13 iterations)
## Equivalent Degrees of Freedom (Df): 45.4384
## Penalized Criterion (RSS): 3842872145
## GCV: 25780520
# Pronostico
predict(spl,232:239)
## $x
## [1] 232 233 234 235 236 237 238 239
## 
## $y
## [1]  5124.1672  3337.0701  1549.9731  -237.1239 -2024.2209 -3811.3179 -5598.4149
## [8] -7385.5119

El ancho de ventana que propone la función por GCV es mejor que el spar de 0,1 con el que hemos realizado el primer spline.

  1. las series de fourier. Justifique el numero de armónicos elegido.
# Periodograma
library(descomponer)
periodograma(jubilados_d)
##          omega frecuencia   periodos    densidad
## 2   0.02719994          1 231.000000  6170135.49
## 3   0.05439987          2 115.500000  4958128.54
## 4   0.08159981          3  77.000000 13233728.59
## 5   0.10879975          4  57.750000 26078015.92
## 6   0.13599968          5  46.200000   575142.89
## 7   0.16319962          6  38.500000   958222.30
## 8   0.19039955          7  33.000000  5808036.93
## 9   0.21759949          8  28.875000  5540132.12
## 10  0.24479943          9  25.666667 10087433.14
## 11  0.27199936         10  23.100000  1400505.56
## 12  0.29919930         11  21.000000  5473326.25
## 13  0.32639924         12  19.250000  1187519.30
## 14  0.35359917         13  17.769231   109952.02
## 15  0.38079911         14  16.500000  1576774.23
## 16  0.40799905         15  15.400000  1535520.04
## 17  0.43519898         16  14.437500  1104286.74
## 18  0.46239892         17  13.588235  3290367.05
## 19  0.48959886         18  12.833333 12803201.00
## 20  0.51679879         19  12.157895 84194821.93
## 21  0.54399873         20  11.550000  2710117.59
## 22  0.57119866         21  11.000000  2951765.31
## 23  0.59839860         22  10.500000   953712.42
## 24  0.62559854         23  10.043478  1106337.50
## 25  0.65279847         24   9.625000  5491130.54
## 26  0.67999841         25   9.240000   504567.50
## 27  0.70719835         26   8.884615  2478199.52
## 28  0.73439828         27   8.555556   238598.76
## 29  0.76159822         28   8.250000   462868.42
## 30  0.78879816         29   7.965517  2818653.08
## 31  0.81599809         30   7.700000   677677.02
## 32  0.84319803         31   7.451613   155188.88
## 33  0.87039796         32   7.218750   760641.88
## 34  0.89759790         33   7.000000  3992155.07
## 35  0.92479784         34   6.794118 11188384.65
## 36  0.95199777         35   6.600000   452327.47
## 37  0.97919771         36   6.416667  1483615.60
## 38  1.00639765         37   6.243243  3812506.92
## 39  1.03359758         38   6.078947 29036658.37
## 40  1.06079752         39   5.923077  2165441.09
## 41  1.08799746         40   5.775000  1206787.76
## 42  1.11519739         41   5.634146  5301287.19
## 43  1.14239733         42   5.500000  1514792.63
## 44  1.16959726         43   5.372093   856850.24
## 45  1.19679720         44   5.250000   199789.62
## 46  1.22399714         45   5.133333  1750394.09
## 47  1.25119707         46   5.021739   543618.11
## 48  1.27839701         47   4.914894  1367344.68
## 49  1.30559695         48   4.812500  2946877.14
## 50  1.33279688         49   4.714286  1322059.80
## 51  1.35999682         50   4.620000   108235.56
## 52  1.38719676         51   4.529412  1492147.27
## 53  1.41439669         52   4.442308  1506692.41
## 54  1.44159663         53   4.358491  1986254.67
## 55  1.46879657         54   4.277778   464231.50
## 56  1.49599650         55   4.200000   336161.95
## 57  1.52319644         56   4.125000  5007647.46
## 58  1.55039637         57   4.052632  1718359.30
## 59  1.57759631         58   3.982759 49168162.07
## 60  1.60479625         59   3.915254   123433.85
## 61  1.63199618         60   3.850000  2747471.41
## 62  1.65919612         61   3.786885  4256029.11
## 63  1.68639606         62   3.725806  2695627.13
## 64  1.71359599         63   3.666667  1199373.57
## 65  1.74079593         64   3.609375  1115042.49
## 66  1.76799587         65   3.553846   190712.79
## 67  1.79519580         66   3.500000  2210532.23
## 68  1.82239574         67   3.447761    36846.25
## 69  1.84959567         68   3.397059   934835.94
## 70  1.87679561         69   3.347826   861704.47
## 71  1.90399555         70   3.300000   507747.74
## 72  1.93119548         71   3.253521  3243362.58
## 73  1.95839542         72   3.208333  3455257.78
## 74  1.98559536         73   3.164384  3026533.41
## 75  2.01279529         74   3.121622  1377235.70
## 76  2.03999523         75   3.080000   283241.29
## 77  2.06719517         76   3.039474   931466.78
## 78  2.09439510         77   3.000000 23572653.87
## 79  2.12159504         78   2.961538   178640.33
## 80  2.14879498         79   2.924051  2783559.44
## 81  2.17599491         80   2.887500   405636.08
## 82  2.20319485         81   2.851852   946867.15
## 83  2.23039478         82   2.817073  1197202.20
## 84  2.25759472         83   2.783133   494184.62
## 85  2.28479466         84   2.750000   851198.91
## 86  2.31199459         85   2.717647   503680.43
## 87  2.33919453         86   2.686047  1879991.19
## 88  2.36639447         87   2.655172  3177665.75
## 89  2.39359440         88   2.625000   158182.25
## 90  2.42079434         89   2.595506    18911.17
## 91  2.44799428         90   2.566667   516338.63
## 92  2.47519421         91   2.538462    42593.69
## 93  2.50239415         92   2.510870  3416666.18
## 94  2.52959408         93   2.483871  1812817.63
## 95  2.55679402         94   2.457447  8475143.91
## 96  2.58399396         95   2.431579  3739893.91
## 97  2.61119389         96   2.406250 29251682.24
## 98  2.63839383         97   2.381443  2125502.59
## 99  2.66559377         98   2.357143  7749052.95
## 100 2.69279370         99   2.333333   347757.54
## 101 2.71999364        100   2.310000  1047964.47
## 102 2.74719358        101   2.287129   175860.44
## 103 2.77439351        102   2.264706  2485511.46
## 104 2.80159345        103   2.242718   430509.34
## 105 2.82879339        104   2.221154  3298473.64
## 106 2.85599332        105   2.200000   819847.40
## 107 2.88319326        106   2.179245   524705.01
## 108 2.91039319        107   2.158879  1587772.68
## 109 2.93759313        108   2.138889   679141.99
## 110 2.96479307        109   2.119266  5755623.84
## 111 2.99199300        110   2.100000   396218.37
## 112 3.01919294        111   2.081081   935679.77
## 113 3.04639288        112   2.062500  4435821.71
## 114 3.07359281        113   2.044248  1970352.87
## 115 3.10079275        114   2.026316 12075325.33
## 116 3.12799269        115   2.008696  8433594.24
gperiodograma(jubilados_d)

# Periodograma acumulado
cpgram(jubilados_d)

td(jubilados_d)
##             s2          min       max
## 1   0.01218869 -0.116959216 0.1343505
## 2   0.02198314 -0.108263564 0.1430462
## 3   0.04812548 -0.099567912 0.1517418
## 4   0.09964085 -0.090872260 0.1604375
## 5   0.10077701 -0.082176608 0.1691331
## 6   0.10266991 -0.073480955 0.1778288
## 7   0.11414330 -0.064785303 0.1865244
## 8   0.12508746 -0.056089651 0.1952201
## 9   0.14501450 -0.047393999 0.2039157
## 10  0.14778111 -0.038698347 0.2126114
## 11  0.15859330 -0.030002695 0.2213070
## 12  0.16093916 -0.021307042 0.2300027
## 13  0.16115636 -0.012611390 0.2386983
## 14  0.16427117 -0.003915738 0.2473940
## 15  0.16730449  0.004779914 0.2560897
## 16  0.16948594  0.013475566 0.2647853
## 17  0.17598583  0.022171218 0.2734810
## 18  0.20127770  0.030866871 0.2821766
## 19  0.36759891  0.039562523 0.2908723
## 20  0.37295257  0.048258175 0.2995679
## 21  0.37878358  0.056953827 0.3082636
## 22  0.38066758  0.065649479 0.3169592
## 23  0.38285307  0.074345132 0.3256549
## 24  0.39370043  0.083040784 0.3343505
## 25  0.39469717  0.091736436 0.3430462
## 26  0.39959269  0.100432088 0.3517418
## 27  0.40006403  0.109127740 0.3604375
## 28  0.40097839  0.117823392 0.3691331
## 29  0.40654645  0.126519045 0.3778288
## 30  0.40788516  0.135214697 0.3865244
## 31  0.40819172  0.143910349 0.3952201
## 32  0.40969432  0.152606001 0.4039157
## 33  0.41758055  0.161301653 0.4126114
## 34  0.43968246  0.169997305 0.4213070
## 35  0.44057600  0.178692958 0.4300027
## 36  0.44350678  0.187388610 0.4386983
## 37  0.45103813  0.196084262 0.4473940
## 38  0.50839810  0.204779914 0.4560897
## 39  0.51267578  0.213475566 0.4647853
## 40  0.51505971  0.222171218 0.4734810
## 41  0.52553205  0.230866871 0.4821766
## 42  0.52852442  0.239562523 0.4908723
## 43  0.53021707  0.248258175 0.4995679
## 44  0.53061174  0.256953827 0.5082636
## 45  0.53406953  0.265649479 0.5169592
## 46  0.53514341  0.274345132 0.5256549
## 47  0.53784451  0.283040784 0.5343505
## 48  0.54366586  0.291736436 0.5430462
## 49  0.54627750  0.300432088 0.5517418
## 50  0.54649132  0.309127740 0.5604375
## 51  0.54943895  0.317823392 0.5691331
## 52  0.55241532  0.326519045 0.5778288
## 53  0.55633904  0.335214697 0.5865244
## 54  0.55725609  0.343910349 0.5952201
## 55  0.55792016  0.352606001 0.6039157
## 56  0.56781243  0.361301653 0.6126114
## 57  0.57120693  0.369997305 0.6213070
## 58  0.66833534  0.378692958 0.6300027
## 59  0.66857917  0.387388610 0.6386983
## 60  0.67400662  0.396084262 0.6473940
## 61  0.68241412  0.404779914 0.6560897
## 62  0.68773915  0.413475566 0.6647853
## 63  0.69010843  0.422171218 0.6734810
## 64  0.69231112  0.430866871 0.6821766
## 65  0.69268786  0.439562523 0.6908723
## 66  0.69705462  0.448258175 0.6995679
## 67  0.69712741  0.456953827 0.7082636
## 68  0.69897411  0.465649479 0.7169592
## 69  0.70067635  0.474345132 0.7256549
## 70  0.70167937  0.483040784 0.7343505
## 71  0.70808642  0.491736436 0.7430462
## 72  0.71491205  0.500432088 0.7517418
## 73  0.72089076  0.509127740 0.7604375
## 74  0.72361140  0.517823392 0.7691331
## 75  0.72417092  0.526519045 0.7778288
## 76  0.72601097  0.535214697 0.7865244
## 77  0.77257717  0.543910349 0.7952201
## 78  0.77293006  0.552606001 0.8039157
## 79  0.77842879  0.561301653 0.8126114
## 80  0.77923010  0.569997305 0.8213070
## 81  0.78110057  0.578692958 0.8300027
## 82  0.78346556  0.587388610 0.8386983
## 83  0.78444179  0.596084262 0.8473940
## 84  0.78612328  0.604779914 0.8560897
## 85  0.78711827  0.613475566 0.8647853
## 86  0.79083206  0.622171218 0.8734810
## 87  0.79710933  0.630866871 0.8821766
## 88  0.79742181  0.639562523 0.8908723
## 89  0.79745916  0.648258175 0.8995679
## 90  0.79847916  0.656953827 0.9082636
## 91  0.79856330  0.665649479 0.9169592
## 92  0.80531269  0.674345132 0.9256549
## 93  0.80889379  0.683040784 0.9343505
## 94  0.82563587  0.691736436 0.9430462
## 95  0.83302378  0.700432088 0.9517418
## 96  0.89080851  0.709127740 0.9604375
## 97  0.89500730  0.717823392 0.9691331
## 98  0.91031503  0.726519045 0.9778288
## 99  0.91100200  0.735214697 0.9865244
## 100 0.91307219  0.743910349 0.9952201
## 101 0.91341959  0.752606001 1.0039157
## 102 0.91832955  0.761301653 1.0126114
## 103 0.91917999  0.769997305 1.0213070
## 104 0.92569590  0.778692958 1.0300027
## 105 0.92731546  0.787388610 1.0386983
## 106 0.92835198  0.796084262 1.0473940
## 107 0.93148852  0.804779914 1.0560897
## 108 0.93283012  0.813475566 1.0647853
## 109 0.94419996  0.822171218 1.0734810
## 110 0.94498267  0.830866871 1.0821766
## 111 0.94683104  0.839562523 1.0908723
## 112 0.95559371  0.848258175 1.0995679
## 113 0.95948601  0.856953827 1.1082636
## 114 0.98334000  0.865649479 1.1169592
## 115 1.00000000  0.874345132 1.1256549
gtd(jubilados_d)

El periodograma apunta a numeros ciclos de alta y baja frecuencia, pero si lo que queremos es tener una señal suavizada, lo mas conveniente es considerar el ciclo de mas baja frecuencía que correspondería al armónico 20.

n=length(jubilados_d)
library(fda)
## Loading required package: splines
## Loading required package: fds
## Loading required package: rainbow
## Loading required package: MASS
## Loading required package: pcaPP
## Loading required package: RCurl
## Loading required package: deSolve
## 
## Attaching package: 'fda'
## The following object is masked from 'package:forecast':
## 
##     fourier
## The following object is masked from 'package:graphics':
## 
##     matplot
fourier = create.fourier.basis(rangeval = c(1,n), nbasis = 20)
plot(fourier)

fourier.fd = smooth.basis( argvals=1:n, y = c(jubilados_d), fdParobj = fourier)
plot(ts(jubilados_d), main="Jubiaciones mensuales")
lines(fourier.fd,col="red")

# Matriz de Regresores
x=t(fourier.fd$y2cMap)
head(x)
##           const        sin1       cos1        sin2       cos2        sin3
## [1,] 0.06042132 0.002334013 0.08541677 0.004666284 0.08532114 0.006995072
## [2,] 0.06049649 0.004882134 0.08541868 0.009749227 0.08501026 0.014586284
## [3,] 0.06071833 0.007432121 0.08555846 0.014813148 0.08462986 0.022092332
## [4,] 0.06107605 0.009981660 0.08582096 0.019842151 0.08416676 0.029461753
## [5,] 0.06155233 0.012528190 0.08618191 0.024819974 0.08359941 0.036643351
## [6,] 0.06212437 0.015068944 0.08660930 0.029730108 0.08289931 0.043586648
##            cos3        sin4       cos4       sin5       cos5       sin6
## [1,] 0.08516185 0.009318641 0.08493900 0.01163526 0.08465277 0.01394319
## [2,] 0.08433094 0.019378400 0.08338278 0.02411080 0.08216865 0.02876891
## [3,] 0.08308918 0.029219610 0.08094685 0.03614594 0.07821736 0.04282363
## [4,] 0.08143179 0.038723622 0.07764881 0.04751518 0.07286315 0.05572948
## [5,] 0.07934887 0.047775054 0.07350967 0.05800470 0.06619088 0.06713873
## [6,] 0.07682677 0.056263293 0.06855473 0.06741641 0.05830529 0.07674322
##            cos6       sin7       cos7       sin8       cos8       sin9
## [1,] 0.08430337 0.01624072 0.08389106 0.01852612 0.08341615 0.02079771
## [2,] 0.08069223 0.03333834 0.07895801 0.03780501 0.07697123 0.04215512
## [3,] 0.07491917 0.04920667 0.07107461 0.05525102 0.06670971 0.06091494
## [4,] 0.06713217 0.06326646 0.06052452 0.07003413 0.05311941 0.07594973
## [5,] 0.05752920 0.07500399 0.04768643 0.08145099 0.03684647 0.08635667
## [6,] 0.04635367 0.08398990 0.03302082 0.08895842 0.01866482 0.09151185
##             cos9      sin10       cos10
## [1,] 0.082878989 0.02305377  0.08227998
## [2,] 0.074737924 0.04637523  0.07226487
## [3,] 0.061853997 0.06615922  0.05654035
## [4,] 0.045005601 0.08094063  0.03628035
## [5,] 0.025211862 0.08962669  0.01300003
## [6,] 0.003671281 0.09157798 -0.01155699

Red Neuronal.

  1. Entrene un arbol de regresión para simular la serie de jubilados.

Se utiliza la función enbed para entrenar el arbol, utilizando 3 meses de desfase, al tratarse de una serie con poca volatilidad.

library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
lag_order <- 3 # número de meses de desfase 
horizon <- 8 # horizonte de predicción (8 meses)

x1 <- embed(jubilados, lag_order + 1) # embedding magic!

y_train <- x1[, 1] # the target
X_train <- x1[, -1]


fit_rf <- randomForest(X_train, y_train)
  
  # Realizo los pronósticos para la muestra de test
  forecasts_rf<- predict(fit_rf, X_train)
 # Graficos
plot(ts(y_train))
lines(ts(c(forecasts_rf)),col="red")

  1. Entrene una red neuronal utilizando neuralnet para simular la serie de jubilados

Se utiliza la función enbed para entrenar la red, utilizando 3 meses de desfase, al tratarse de una serie con poca volatilidad.

y_train <- x1[, 1] # the target
X_train <- x1[, -1]
datos_rn=data.frame(y_train,X_train)
# NORMALIZACION DE VARIABLES
datos_nrm <- scale(datos_rn)
library(neuralnet)
## 
## Attaching package: 'neuralnet'
## The following object is masked from 'package:RJDemetra':
## 
##     compute
mod3 <- neuralnet(y_train ~., data=datos_nrm)
#plot(mod3,rep="best")
test <- data.frame(datos_nrm[,-1])
myprediction<- predict(mod3, test)
plot(ts(datos_nrm[,1]))
lines(ts(myprediction),col="red")

  1. Estima una red neuronal con grnn_forecasting. Realice un pronóstico para lo que queda de año.

Se entrena la red con solo tres desfases, como en los casos anteriores.

library(tsfgrnn)
pred <- grnn_forecasting(c(jubilados),h = 8, lags = 1:3)
pred$prediction
## Time Series:
## Start = 233 
## End = 240 
## Frequency = 1 
## [1] 6356081 6363594 6371176 6378741 6386310 6393878 6401446 6409014
plot(pred)

  1. Utilice la función knn_forecasting, para pronosticar .

Se utilizan 3 desfases al igual que en las redes anteriores.

library(tsfknn)
## 
## Attaching package: 'tsfknn'
## The following object is masked from 'package:tsfgrnn':
## 
##     rolling_origin
pred <- knn_forecasting(c(jubilados), h = 12, lags = 1:3,k=3)
plot(pred)