Actividad 2
Ejercicio 1: serie Tasa Actividad
datos1=read.csv2("Tasa actividad, paro y empleo Alicante.csv")
#------------SEGMENTACIÓN DE DATOS-------------------------------------
datos1_tasa.actividad= datos1 %>%
filter(Tasas=="Tasa de actividad")
datos1_tasa.empleo= datos1 %>%
filter(Tasas=="Tasa de empleo de la poblaci\xf3n")
datos1_tasa.paro= datos1 %>%
filter(Tasas=="Tasa de paro de la poblaci\xf3n")
#----------------------------------------------------------------------
transpose1 <- t(datos1_tasa.actividad)
transpose1 <- as.data.frame(transpose1)
rev_data_frame1 <- rev(transpose1)
rev_data_frame1 <- t(rev_data_frame1)
rev_data_frame1 <- as.data.frame(rev_data_frame1)
datos1_tasa.actividad=rev_data_frame1
transpose2 <- t(datos1_tasa.empleo)
transpose2 <- as.data.frame(transpose2)
rev_data_frame2 <- rev(transpose2)
rev_data_frame2 <- t(rev_data_frame2)
rev_data_frame2 <- as.data.frame(rev_data_frame2)
datos1_tasa.empleo=rev_data_frame2
transpose3<- t(datos1_tasa.paro)
transpose3 <- as.data.frame(transpose3)
rev_data_frame3 <- rev(transpose3)
rev_data_frame3 <- t(rev_data_frame3)
rev_data_frame3 <- as.data.frame(rev_data_frame3)
datos1_tasa.paro=rev_data_frame3Tasa Actividad
knitr::kable(head(datos1_tasa.actividad))| Sexo | Provincias | Tasas | Periodo | Total | |
|---|---|---|---|---|---|
| V84 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2002T1 | 54.40 |
| V83 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2002T2 | 55.51 |
| V82 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2002T3 | 55.17 |
| V81 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2002T4 | 55.19 |
| V80 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2003T1 | 55.11 |
| V79 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2003T2 | 54.98 |
datos2_tasa.paro=datos1_tasa.paro
n=nrow(datos1_tasa.paro)
Tasas=rep(c("Tasa de paro"),each=n)
datos2_tasa.paro$Tasas=Tasas
comprobacion=datos2_tasa.paro==datos1_tasa.paro
datos1_tasa.paro=datos2_tasa.paro
#knitr::kable(head(datos1_tasa.paro))A) Análisis básico
filas=nrow(datos1_tasa.actividad)
ultimovalor=datos1_tasa.actividad[filas,];ultimovalor| Sexo | Provincias | Tasas | Periodo | Total | |
|---|---|---|---|---|---|
| V1 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2022T4 | 60.21 |
- Las serie termina el cuarto trimestre de 2022
filter(datos1_tasa.actividad,Periodo=="2014T1")| Sexo | Provincias | Tasas | Periodo | Total | |
|---|---|---|---|---|---|
| V36 | Ambos sexos | 03 Alicante/Alacant | Tasa de actividad | 2014T1 | 57.49 |
- El primer cuatrimestre de 2014 tenemos un valor de 57.49
B) Visualización básica de la serie
- En este gráfico observamos una visualización básica de la gráfica donde de primera mano podemos detectar algunos puntos importantes
actividad=ts(as.numeric(datos1_tasa.actividad[,"Total"]),
start=c(2002,1),
frequency=4)
autoplot(actividad,
xlab="",
ylab="Actividad",
main="Actividad por años")- En este gráfico tratamos de observar la tendencia y los ciclos de la serie de Actividad, hemos marcado los puntos más relevantes de la serie para observar sucesos importantes (por ejemplo el Covid en 2020)
actividadAnual=aggregate(actividad,FUN=sum)
min=min(actividadAnual)
min=max(actividadAnual)
autoplot(actividadAnual,
xlab="",
ylab="Actividad",
main="Actividad por años") +
geom_hline(yintercept = c(220.27,225.02,240.09),linetype="dashed",color="red")+
geom_vline(xintercept=c(2002,2009,2020),linetype="dashed",color="red")- En este gráfico, observamos como no hay ninguna tendencia clara, pues los puntos del gráfico se reparten más bien de forma aleatoria por todo el “mapa”
actividadDesv=aggregate(actividad,FUN=sd)
plot(as.numeric(actividadAnual),as.numeric(actividadDesv))- En este gráfico comparamos la Actividad por trimestres y años; de forma clara, vemos que los últimos trimestres del año presentan de media más actividad que los primeros, además, podemos ver que el año anterior fue un año de baja actividad y que este está siendo una recuperación, particularmente para el segundo cuatrimestre, donde en el año anterior se obtuvieron datos catastróficos.
actividadb=window(actividad,start=2002)
ggsubseriesplot(actividadb)+
ylab("Actividad")+
xlab("")+
ggtitle("")- Finalmente en este gráfico pretendemos buscar observaciones anómalas en los años comparando la actividad total entre todos ellos. Aún viendo que no hay ninguna, el gráfico nos ayuda a ver que la Activida tiende a crecer, pues los registros más altos son 2022, 2008,2009 y 2019, no estando 2020 y 2021 por el Covid.
ggseasonplot(actividadb,
year.labels=TRUE,
xlab="",
ylab="Actividad",
main="")C) Regresiones locales
actividad_stl=stl(actividad,
s.window="periodic",
robust=TRUE)
autoplot(actividad_stl,
xlab="",
main="Descomposición de Actividad por regresores locales ponderados")- Observamos en este gráfico como la tendencia actual de la Actividad es mantenerse a niveles medios.
D) Medias móviles
- Gracias a la siguiente gráfica podemos suavizar la tendencia que ha llevado la Actividad
- En 2-Ma estamos tomando la media de dos trimestres, el anterior, y el actual.
- En 4-Ma estamos tomando la media de un año, un trimestre anterior, el actual y dos posteriores.
- En 2*4-Ma estamos tomando la media de dos años, el año anterior y el actual.
actividad_tsibble=tsibble(
time= as.numeric(time(actividad)),
actividad=as.numeric(datos1_tasa.actividad$Total),
index=time)|>
mutate(
`2-Ma` = slider::slide_dbl(actividad,mean, .before=1, .after=0, .complete=FALSE),
`4-Ma` = slider::slide_dbl(actividad,mean,.before=1,.after=2,.complete=FALSE),
`2*4-Ma` = slider::slide_dbl(`4-Ma`,mean,.before=1,.after=0,.complete=FALSE),
)
p=actividad_tsibble |>
ggplot()+
geom_line(aes(x=time,y=actividad),col="gray45",size=0.55)+
geom_line(aes(x=time,y=`2-Ma`),col="red",size=0.65)+
geom_line(aes(x=time,y=`2*4-Ma`),col="green",size=0.65)+
xlab("año")+
scale_color_manual(values=c("grey45","red","green"),labels=c("A","B","C"))
p + theme(legend.position="left")E) Predicción de la Serie Actividad
- Predicción de la serie bajo el método Ingenuo II en los próximos dos años
derivaActividad=rwf(actividad,h=8,drift=TRUE)
knitr::kable(derivaActividad)| Point Forecast | Lo 80 | Hi 80 | Lo 95 | Hi 95 | |
|---|---|---|---|---|---|
| 2023 Q1 | 60.28 | 58.95313 | 61.60687 | 58.25073 | 62.30927 |
| 2023 Q2 | 60.35 | 58.46239 | 62.23761 | 57.46315 | 63.23685 |
| 2023 Q3 | 60.42 | 58.09460 | 62.74540 | 56.86361 | 63.97639 |
| 2023 Q4 | 60.49 | 57.78929 | 63.19071 | 56.35962 | 64.62038 |
| 2024 Q1 | 60.56 | 57.52321 | 63.59679 | 55.91563 | 65.20437 |
| 2024 Q2 | 60.63 | 57.28452 | 63.97548 | 55.51353 | 65.74647 |
| 2024 Q3 | 60.70 | 57.06622 | 64.33378 | 55.14262 | 66.25738 |
| 2024 Q4 | 60.77 | 56.86381 | 64.67619 | 54.79599 | 66.74401 |
autoplot(actividad,
series="Actividad",
xlab="",
ylab="Actividad",
main="Predicción sin estacionalidad")+
autolayer(derivaActividad,series="Ingenuo II", PI=FALSE)+
labs(colour="Métodos")+
theme(legend.position=c(0.1,0.8))snaive.Actividad=snaive(actividad,h=8,level=95)
autoplot(snaive.Actividad,
xlab="",
ylab="Actividad",
main="Predicción con estacionalidad",
PI=FALSE)- Utilizamos el método sin estacionalidad porque cuando lo analizamos vemos que este ofrece mejor error porcentual (MAPE=1.355733 frente al estacional MAPE=2.094922), además, el intervalo de confianza de las predicciones es mucho más fiable (ACF1=-0.1053811 frente a ACF1=0.5038194)
datosderiva=accuracy(derivaActividad)
datossnaive=accuracy(snaive.Actividad)
comparacion=as.data.frame(rbind(datosderiva,datossnaive))
rownames(comparacion)=c("Sin estacionalidad","Con estacionalidad")
comparacion| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Sin estacionalidad | 0.000000 | 1.022960 | 0.7826506 | -0.0147015 | 1.355733 | 0.6459512 | -0.1053811 |
| Con estacionalidad | 0.163125 | 1.437565 | 1.2116250 | 0.2558152 | 2.094921 | 1.0000000 | 0.5038194 |
Ejercicio 2: Serire Tasa Paro
A) Graficar y comentar la serie y su esquema
paro=ts(as.numeric(datos1_tasa.paro[,"Total"]),
start=c(2002,1),
frequency=4)
autoplot(paro,
xlab="",
ylab="Paro",
main="Paro trimestral")- Podemos observar de una forma clara, como el paro tenía su nivel más bajo en 2005, donde el empleo era pleno y este aumenta en 2007 a raíz de la crisis económica de ese periodo. Vuelve a bajar desde 2013 hasta 2019 cuando por la Covid-19 vuelve a aumentar. Actualmente y desde la pandemia, el paro está bajando.
paroAnual=aggregate(paro,FUN=sum)
min=min(paroAnual)
max=max(paroAnual)
autoplot(paroAnual,
xlab="",
ylab="Paro",
main="Paro anual") +
geom_hline(yintercept = c(max,min),linetype="dashed",color="red")+
geom_vline(xintercept=c(2007,2019),linetype="dashed",color="red")- En este gráfico vemos que se presenta un esquema aditivo por la aleatoriedad de los puntos en el gráfico
paroDesv=aggregate(paro,FUN=sd)
plot(as.numeric(paroAnual),as.numeric(paroDesv))B) Predicción de ingenuo con estacionalidad para los próximos dos años
snaive.Paro=snaive(paro,h=8,level=95)
autoplot(snaive.Paro,
xlab="",
ylab="Paro",
main="Predicción con estacionalidad",
PI=FALSE)knitr::kable(snaive.Paro)| Point Forecast | Lo 95 | Hi 95 | |
|---|---|---|---|
| 2023 Q1 | 14.78 | 8.603176 | 20.95682 |
| 2023 Q2 | 14.13 | 7.953176 | 20.30682 |
| 2023 Q3 | 13.61 | 7.433176 | 19.78682 |
| 2023 Q4 | 15.21 | 9.033176 | 21.38682 |
| 2024 Q1 | 14.78 | 6.044652 | 23.51535 |
| 2024 Q2 | 14.13 | 5.394652 | 22.86535 |
| 2024 Q3 | 13.61 | 4.874652 | 22.34535 |
| 2024 Q4 | 15.21 | 6.474652 | 23.94535 |
C) Aplicación del método Holt-Winters
- Ofrecemos la vista de ambas posibilidades para mostrar la similitud entre ambos modelos, pero nos debemos quedar con el esquema aditivo, pues en la graficación de puntos no observamos ninguna tendencia creciente.
# Aditivo
paropred <- window(paro, start = c(2002,1), end = c(2022, 4))
paroets <- ets(paropred,
model = "AAA",
damped = FALSE)
aditivo=data.frame(alpha=paroets$par["alpha"],beta=paroets$par["beta"],gamma=paroets$par["gamma"],row.names=c("Aditivo"))
# Multiplicativo
parob <- window(paro, start = 2002)
parobEts <- ets(parob,
model = "MAM",
damped = FALSE)
multiplicativo=data.frame(alpha=parobEts$par["alpha"],beta=parobEts$par["beta"],gamma=parobEts$par["gamma"],row.names=c("Multiplicativo"))
rbind(aditivo,multiplicativo)| alpha | beta | gamma | |
|---|---|---|---|
| Aditivo | 0.9789324 | 0.1034869 | 0.0001000 |
| Multiplicativo | 0.8939725 | 0.1469861 | 0.0001001 |
- Los valores bajos en Beta (0.103) y Gamma (0.0001) nos indica que la pendiente y la estacionalidad se modifican muy lentamente. Sin embargo, el Alfa tan grande (0.9789324) indica que la serie varía constantemente en el tiempo
D) Analizar los errores del modelo Aditivo
- Identificamos periodos cruciales como puede ser 2009, donde el paro ascendió de forma desorbitada, al igual que en el primer trimestre de 2020, con el coronavirus.
error=residuals(paroets)
sderror=sd(error)
autoplot(error,series="Error",
colour="black",
xlab="Semana",
ylab="Error",
main="")+
geom_hline(yintercept= c(-3,-2,2,3)*sderror,
colour=c("red","blue","blue","red"),lty=2)+
scale_x_continuous(breaks=seq(6,26,2))+
geom_vline(xintercept=c(2009,2020.25),linetype="dashed",color="red")E) Predicción método Holt-Winters aditivo
parodf=forecast(paroets,
h=8,
level=95)
knitr::kable(parodf)| Point Forecast | Lo 95 | Hi 95 | |
|---|---|---|---|
| 2023 Q1 | 15.26319 | 12.247755 | 18.27862 |
| 2023 Q2 | 14.42260 | 9.978919 | 18.86628 |
| 2023 Q3 | 13.59436 | 7.890484 | 19.29824 |
| 2023 Q4 | 14.29876 | 7.395760 | 21.20176 |
| 2024 Q1 | 14.37313 | 6.292589 | 22.45366 |
| 2024 Q2 | 13.53254 | 4.277539 | 22.78754 |
| 2024 Q3 | 12.70430 | 2.267422 | 23.14118 |
| 2024 Q4 | 13.40870 | 1.776409 | 25.04098 |
autoplot(paro,
xlab="",
ylab="Paro",
main="Predicción paro Holt-Winters aditivo",
colour="black",
PI=FALSE)+
autolayer(parodf,series="Holt-Winters aditivo",PI=FALSE)+
autolayer(snaive.Paro,series="Ingenuo con Estacionalidad", PI=FALSE)+
labs(colour="Métodos")+
theme(legend.position=c(0.1,0.8))- Observamos que la serie predecida con Ingenuo con tendencia no tiene un intervalo de confianza de las predicciones fiable (ACF1=0.804918751). Además, tiene un error porcentual menor la predicción Holt-Winters (MAPE=6.702592%).
datosHolt=accuracy(parodf)
datossnaiveParo=accuracy(snaive.Paro)
comparacion=as.data.frame(rbind(datosHolt,datossnaiveParo))
rownames(comparacion)=c("Holt-Winters Aditivo","Ingenuo con Estacionalidad")
comparacion| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Holt-Winters Aditivo | -0.1193827 | 1.463420 | 1.097746 | -0.8978805 | 6.702592 | 0.4634041 | 0.0053174 |
| Ingenuo con Estacionalidad | 0.1873750 | 3.151499 | 2.368875 | -0.1132280 | 13.362802 | 1.0000000 | 0.8049188 |
F) Validación Training Set
- Hemos realizado la predicción del modelo para Holt-Winters Aditivo y para Ingenuo con Estacionalidad
- Aunque ambos tengan un error alto, el que mejor se adecua es Holt-Winters Aditivo (MAPE=6.595903%),
#Eliminamos los últimos 14 trimestres
paroIntra=subset(paro,end=length(paro)-14)
# additive predicción para esos 14 trimestres para Holt-Winters aditivo
paroIntraEts=ets(paroIntra,
mode="AAA",
damped=FALSE)
paroEtsExtraPre=forecast(paroIntraEts,
h=14,
level=95)
# Hacemos predicción para Ingenuo con Estacionalidad
snaive.paroIntra=snaive(paroIntra,h=14,level=95)
#Información
informacionmodelos=data.frame(rbind(accuracy(snaive.paroIntra),accuracy(paroEtsExtraPre)))
rownames(informacionmodelos)=c("Ingenuo II Estacionalidad","Holt-Winters Aditivo")
informacionmodelos| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Ingenuo II Estacionalidad | 0.2754545 | 3.125400 | 2.297273 | 0.6106416 | 12.693378 | 1.0000000 | 0.8191554 |
| Holt-Winters Aditivo | -0.0386868 | 1.373433 | 1.067346 | 0.0436654 | 6.595903 | 0.4646142 | 0.0871120 |
#Gráfica
autoplot(paro, series="Paro",linetype="solid")+
autolayer(paroIntra,series="Reducida",linetype="dashed")+
autolayer(paroEtsExtraPre,series="Holt-Winter Aditivo",linetype="solid",PI=F)+
autolayer(snaive.paroIntra,series="Ingenuo Estacionalidad",linetype="solid",PI=F)Ejercicio 3: Serie Consumo de Alimentos por Cápita
A) Análisis Básico
- La serie termina en el año 2021
- Llama mucho la atención 2019 y 2020, donde por el Coronavirus se produce este movimiento tan extraño en la serie, donde pasamos de 606.64 a 671.82 en tan solo un año
- Como se comenta en el apartado anterior, la explicación es la Covid-19
consumo=read.csv2("alimentacionpc.csv")
consumo=ts(consumo,
start=c(1987),
frequency=1)B) Graficar las componentes de la serie
- En este gráfico vemos los cortes que presenta el gráfico en sus momentos más cruciales (1993 - 1995) una bajada drástica en el Consumo de Alimentos por Cápita y (2019 - 2020) una drástica subida
autoplot(consumo,
xlab="",
ylab="Consumo de Alimentos por Cápita",
main="Consumo anual")+
geom_vline(xintercept=c(2019,2020,1993,1995),linetype="dashed",color="red")- Con este gráfico vemos como claramente estamos ante un esquema multiplicativo, pues los puntos dispersos en el gráfico muestran una tendencia creciente
- Al tratarse de una serie anual, no podremos estudiar la estacionalidad de la serie
consumoDesv=aggregate(consumo,FUN=sd)
plot(as.numeric(consumo),as.numeric(consumoDesv))C) Estacionariedad de la Serie
- La estacionariedad hace referencia a la característica de una serie que la hace estacionaria en el tiempo, es decir, que el transcurso del tiempo no influye en la variación de los valores de la misma.
- Según la definición anterior podemos decir que la serie que estamos trabajando no es estacionaria, pues podemos ver como el paso del tiempo si influye en el aumento o la reducción de los valores de consumo.
D) Predicción de valores con método de la media
mediaConsumo=meanf(consumo,h=2)
autoplot(consumo,
seires="Consumo Anual per Cápita",
xlab="",
ylab="Consumo",
main="Consumo Anual de Alimentos per Cápita")+
autolayer(mediaConsumo,series="Media",PI=FALSE)+
labs(colour="Métodos")+
theme(legend.position=c(0,0.8))- Podemos remarcar que las predicciones no presentan un intervalo de confianza muy fiable (ACF1=0.292032).
info=as.data.frame(accuracy(mediaConsumo))
rownames(info)="Método de la media"
info| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Método de la media | 0 | 16.33048 | 13.07631 | -0.0649985 | 2.04427 | 1.068095 | 0.292032 |
E) Alisado Exponencial Simple
Vemos que alpha es igual a 0.3742, por lo que la serie no varía muy drásticamente con el tiempo, además, el error porcentual (MAPE=1.930308) es realmente bajo indicando así un ajuste bueno
Además, \(y_{t+2}=y_{t+1}=l_t=633.055\)
\(l_{t}=\alpha*y_t+(1-\alpha)l_{t-1}\)
consumo.alisado=window(consumo,start=c(1987),end=c(2021))
consumoEts=ets(consumo.alisado,model="ANN")
alpha=consumoEts$par["alpha"]
lt=consumoEts$states[36,]
info2=data.frame(Alpha=alpha,Lt=lt,"Lt-1"=653.9005)
rownames(info2)="Datos"
info2| Alpha | Lt | Lt.1 | |
|---|---|---|---|
| Datos | 0.3741634 | 633.055 | 653.9005 |
F) Alisado Exponencial Simple
- Alpha es un valor entre 0 y 1, significando lo siguiente para cada caso:
- Si Alpha=0, la serie permanece constante en el tiempo
- Si Alpha=1, la serie varía constantemente en el tiempo
- En nuestro caso, como se ha comentado antes, alpha=0.3742, por lo que podemos decir que se mantiene medianamente constante en el tiempo.
consumodf=forecast(consumoEts,h=2,level=95)
autoplot(consumodf,xlab="",ylab="Consumo per cápita",main="Predicción del consumo per cápita")G) Comparación de ambas predicciones
- Gráficamente podemos ver como ambas predicciones son bastante similares
autoplot(consumo,series="Consumo")+
autolayer(mediaConsumo,series="Predicción Media",PI=F)+
autolayer(consumodf,series="Alisado Exponencial Simple",PI=F)- Si vemos los datos más técnicos, vemos como aunque ambos son parecidos en cierto sentido, podemos decir que el Alisado Exponencial Simple tiene un error menor (ACF1=0.0853%), teniendo así un intervalo de confianza aceptable.
infocomparacion=data.frame(rbind(accuracy(mediaConsumo),accuracy(consumodf)))
rownames(infocomparacion)=c("Método Medias","Alisado Exponencial Simple")Ejercicio 4: Serie Población de Niños
Datos
poblacion=read.csv2("Población en miles de personas.csv")
poblacion=poblacion%>%filter(Edad=="De 0 a 4 anos")
poblacion=poblacion%>%filter(Unidad=="Valor absoluto")
poblacion=poblacion%>%filter(Sexo=="Ambos sexos")
poblacion=as.data.frame(apply(poblacion,2,rev))
head(poblacion)| Edad | Sexo | Unidad | Periodo | Total |
|---|---|---|---|---|
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2002T1 | 1945.0 |
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2002T2 | 1963.5 |
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2002T3 | 1982.9 |
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2002T4 | 2007.6 |
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2003T1 | 2027.3 |
| De 0 a 4 anos | Ambos sexos | Valor absoluto | 2003T2 | 2050.5 |
poblacion=ts(as.numeric(poblacion[,"Total"]),
start=c(2002,1),
frequency=4)A) Graficar y comentar análisis básico
Vemos como el punto más alto de niños se encuentra en el último cuatrimestre de 2010, punto a partir del cual la cifra comienza a bajar considerablemente, teniendo un parón en la bajada en el tercer cuatrimestre de 2018.
Este parón en la bajada no es siquiera notable si observamos la serie de forma anual
autoplot(poblacion,
xlab="",
ylab="Poblacion",
main="Poblacion")+
geom_vline(xintercept=c(2010.75,2018.5),linetype="dashed",color="red")+
geom_hline(yintercept=c(as.numeric(max(poblacion)),2084.1),color="blue",linetype="dashed")poblacionAnual=aggregate(poblacion, FUN = sum)
autoplot(poblacionAnual)+
geom_hline(yintercept=c(max(poblacionAnual)),color="red",linetype="dashed")- En el siguiente gráfico observamos el tipo de esquema que sigue la serie, si nos fijamos, los puntos se esparcen de forma “aleatoria” por el gráfico, por lo que es un esquema aditivo.
poblacionDesv <- aggregate(poblacion, FUN = sd)
plot(as.numeric(poblacionAnual ), as.numeric(poblacionDesv))- Como vemos, se trata de una serie que se ha movido exactamente igual a lo largo de cada trimestre en los últimos años, por lo que tiene la misma forma y prácticamente la misma media
ggsubseriesplot(poblacion)+
ylab("Poblacion")+
xlab("Años")+
ggtitle("")- Como vemos en este gráfico y se podía intuir anteriormente, la población en niños no aumenta estacionalmente, sino que se mantiene igual durante todo el año
ggseasonplot(poblacion,
year.labels=TRUE,
xlab="",
ylab="Poblacion",
main="")B) Holt-Winters
- \(\alpha=0.9994\) Lo que indica que esta serie varía mucho con el paso del tiempo.
- \(\beta=0.1913\) indica que la pendiente apenas se mueve con el paso de tiempo, y tiende a mantenerse constante.
- \(\gamma\approx0\) indica que la estacionalidad permanece constante en el tiempo. Esta componente no aporta apenas información
poblacionholt=window(poblacion,start=c(2002,1),end=c(2022,4))
poblacionEts=ets(poblacionholt,
model="AAA",
damped=FALSE)
summary(poblacionEts)## ETS(A,A,A)
##
## Call:
## ets(y = poblacionholt, model = "AAA", damped = FALSE)
##
## Smoothing parameters:
## alpha = 0.9994
## beta = 0.1913
## gamma = 6e-04
##
## Initial states:
## l = 2076.3123
## b = 13.9701
## s = 2.5819 4.433 -3.3095 -3.7054
##
## sigma: 19.5961
##
## AIC AICc BIC
## 881.6368 884.0692 903.5142
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1.823973 18.63957 10.04761 -0.07581853 0.4718947 0.1680589
## ACF1
## Training set 0.01296637
C) Componentes estacionales s1,s2,s3 y s4
- Como terminamos en el último cuatrismetre de 2022, tenemos que s1 es el valor estacional para el último dato (último trimester de 2022), s2 tercer trimestre de 2022, s3 segundo trimestre y s1 primer trimestre.
tail(poblacionEts$states,1)## l b s1 s2 s3 s4
## 2022 Q4 1796.757 -15.34179 2.642959 4.357625 -3.349597 -3.749964
D) Predicción a 5 años
- En este ejercicio he probado a bajar el confidence level para estrechar el canal de predicción y ver otra alternativa. Vemos que este modelo está destinado a seguir bajando, a no ser que se produzca algo independiente a la serie que lo haga ser modificado.
poblaciondf=forecast(poblacionEts,
h=20,
level=70)
autoplot(poblaciondf,
xlab="Población + 5 años predicción",
ylab="Población",
main="Predicción Holt-Winters Aditivo",
PI=TRUE)E) Validación Cruzada
- El error a un trimestre vista es de 7.77, sin embargo, si nos vamos al final del periodo (8 trimestres o 2 años), el error asciende a 58.58
k <- 32 #Minimo numero de datos para estimar
h <- 8 #Horizonte de las prediciciones
TT <- length(poblacion)#Longitud serie
s <- TT - k - h #Total de estimaciones
rmse <- matrix(NA, s + 1, h)
mape <- matrix(NA, s + 1, h)
for (i in 0:s) {
train.set <- subset(poblacion, start = i + 1, end = i + k)
test.set <- subset(poblacion, start = i + k + 1, end = i + k + h)
poblacionEts <- ets(train.set,
model = "AAA",
damped = FALSE)
fcast <- forecast(poblacionEts, h = h)
rmse[i + 1,] <- (test.set - fcast$mean)^2
mape[i + 1,] <- 100*abs(test.set - fcast$mean)/test.set
}
rmse.ets <- sqrt(colMeans(rmse))
mape.ets <- colMeans(mape)
error=data.frame(t(round(rmse.ets, 2)))
colnames(error)=c("1T","2T","3T","4T","5T","6T","7T","8T")
rownames(error)=c("Error")
error| 1T | 2T | 3T | 4T | 5T | 6T | 7T | 8T | |
|---|---|---|---|---|---|---|---|---|
| Error | 7.77 | 14.1 | 20.86 | 26.68 | 32.64 | 41.51 | 49.96 | 58.58 |
Ejercicio 5: Serie Pasajeros
A) Graficar la serie completa y comentar
- Salta claramente a la vista el año 2020, como vemos el año comenzó como cualquier otro, pero en el momento en el que la Covid se hizo una realidad, el transporte bajó drásticamente.
pasajeros=read.csv2("Pasajeros.csv")
pasajeros=ts(pasajeros,
start=c(1996,1),
frequency=12)
pasajerosb=window(pasajeros,start=1996)
ggseasonplot(pasajerosb,
year.labels=T,
xlab="Años",
ylab="Pasajeros",
main="")B) Recortar la serie
- Seguimos con el ejercicio a partir de aquí
pasajeros2=as.data.frame(pasajeros)
pasajeros=pasajeros2[1:(312-24),]
tibble(pasajeros=head(pasajeros))| pasajeros |
|---|
| 216721 |
| 213787 |
| 221582 |
| 204038 |
| 220739 |
| 201742 |
pasajeros=ts(pasajeros,
start=c(1996,1),
frequency=12)C) Graficas y comentar
- En este gráfico apreciamos como el crecimiento de usuarios en el transporte público desde 1996 creción significativamente hasta llegar a un punto (2007) donde se paró el crecimiento y comenzó a bajar. Todo esto hasta 2013 donde de nuevo se retoma el uso de transporte y crece hasta sus niveles más altos en 2019.
pasajerosAnual=aggregate(pasajeros,FUN=sum)
autoplot(pasajerosAnual,
xlab="Años",
ylab="Pasajeros",
main="Pasajeros transporte público 1996-2019")+
geom_vline(xintercept=c(2007,2013,2019),linetype="dashed",color="red")- En el siguiente gráfico observamos el esquema del gráfico, donde vemos que es un esqauema aditivo por la aleatoriedad de los puntos en el gráfico
pasajerosDesv=aggregate(pasajeros,FUN=sd)
pasajerosAnual=aggregate(pasajeros,FUN=sum)
plot(as.numeric(pasajerosAnual),as.numeric(pasajerosDesv))- En este gráfico estudiamos la estacinalidad de los pasajeros en la serie. Como vemos, en los periodos vacacionales el uso de transporte público disminuye drásticamente (de Junio a Septiembre se mantiene bajo, al igual que en Navidad) y se encuentra en sus mayores niveles en momentos de máximo trabajo en el calendario anual (Octubre y Marzo).
pasajerosb=window(pasajeros,start=1996)
ggsubseriesplot(pasajerosb)+
ylab("Pasajeros")+
xlab("")+
ggtitle("Estacionalidad de pasajeros")D) Medias móviles
actividad_tsibble=tsibble(
time= as.numeric(time(pasajeros)),
actividad=as.numeric(pasajeros),
index=time)|>
mutate(
`3-Ma` = slider::slide_dbl(actividad,mean, .before=1, .after=1, .complete=FALSE),
`12-Ma` = slider::slide_dbl(actividad,mean,.before=5,.after=6,.complete=FALSE),
`2*12-Ma` = slider::slide_dbl(`12-Ma`,mean,.before=1,.after=0,.complete=FALSE),
)
p=actividad_tsibble |>
ggplot()+
geom_line(aes(x=time,y=actividad),col="gray45",size=0.55)+
geom_line(aes(x=time,y=`3-Ma`),col="red",size=0.65)+
geom_line(aes(x=time,y=`12-Ma`),col="green",size=0.65)+
geom_line(aes(x=time,y=`2*12-Ma`),col="blue",size=0.65)+
xlab("año")+
scale_color_manual(values=c("grey45","red","green"),labels=c("A","B","C"))
p + theme(legend.position="left")Ejercicio 6: Serie Aforo Oropesa
A) Graficar y comentar las componentes principales
En el siguiente gráfico vemos la evolución anual de los vehículos que circulan por Oropesa, donde observamosdos tramos principales. De 1960 a 1982 donde vemos una montaña de crecimiento en la circulación (se puede deber quizá a una industria muy explotada en la zona que requería mucha circulación de tráfico). Y el siguiente tramo desde 1982 en adelante, donde vemos un crecimiento sustancial de vehículos en la zona. A este aumento podemos atribuirle la expansión del automóvil, pues es evidente que en este tramo el número de vehículos disponibles para la población creción sustancialmente.
Desde 2007 en adelante el número de vehículos en ese km disminuye por la construcción de una autopista que permite una mayor fluidez y velocidad en el tráfico.
vehiculos=read.csv2("aforo_oropesa.csv")
vehiculos=ts(vehiculos,
start=c(1960),
frequency=1)
autoplot(vehiculos,
xlab="Años",
ylab="Vehículos",
main="Vehículos por Oropesa")+
geom_hline(yintercept = c(max(vehiculos)),color="red",linetype="dotdash",lwd=1)+
geom_vline(xintercept=c(1982,2007),linetype="twodash",color="red")- Gracias al siguiente gráfico podemos afirmar que se trata de un esquema multiplicativo debido a la continuidad de los puntos en el gráfico.
vehiculosDesv=aggregate(vehiculos,FUN=sd)
plot(as.numeric(vehiculos),as.numeric(vehiculosDesv))B) Predicción método media y deriva
El método que mejor resultado obtiene es el método Deriva, pero realmente ninguno de los dos ofrece un intervalo de confianza para el error bueno.
La diferencia principal entre ambos métodos es que la media simplemente calcula la media total de la serie para ofrecer la información y el método deriva añade a la observación anterior el incremento medio de la serie * el número de periodos. (Si se quiere obtener predicción de dos años en una serie anual, se calcula el incremento medio y se multiplica por dos)
mediaVehiculos=meanf(vehiculos,h=5)
mediaVehiculosdf=as.data.frame(mediaVehiculos)
# método deriva
derivaVehiculos=rwf(vehiculos, h=5, drift=T)
derivaVehiculosdf=as.data.frame(derivaVehiculos)
predicciones=data.frame(Metodo_Media=mediaVehiculosdf[,"Point Forecast"],Metodo_deriva=derivaVehiculosdf[,"Point Forecast"])
rownames(predicciones)=c(2020:2024)
predicciones| Metodo_Media | Metodo_deriva | |
|---|---|---|
| 2020 | 10379.37 | 9824.746 |
| 2021 | 10379.37 | 9961.492 |
| 2022 | 10379.37 | 10098.237 |
| 2023 | 10379.37 | 10234.983 |
| 2024 | 10379.37 | 10371.729 |
vehiculosderiva=accuracy(derivaVehiculos)
vehiculosmedia=accuracy(mediaVehiculos)
autoplot(vehiculos,
series="Original")+
autolayer(derivaVehiculos,PI=F,series="predicción deriva")+
autolayer(mediaVehiculos,PI=F,series="predicción media")comparacion=as.data.frame(rbind(vehiculosderiva,vehiculosmedia))
rownames(comparacion)=c("Método Deriva","Método Media")
comparacion| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Método Deriva | 0 | 925.7564 | 640.9457 | 0.8655602 | 6.865007 | 0.9444505 | 0.3001172 |
| Método Media | 0 | 3638.2852 | 2826.2300 | -30.5042998 | 50.645009 | 4.1645247 | 0.9188694 |
C) Tipos de Alisado Exponencial posibles
Gracias a las gráficas del apartado “A)”, sabemos que el esquema de la serie es multiplicativo, por lo que solamente podemos aplicar métodos sin estacionalidad (porque es una serie anual) y con error multiplicativo.
Hemos aplicado los métodos MAN (error multiplicativo, tendencia aditiva y sin estacionalidad) y MNN (error multiplicativo, sin tendencia y sin estacionalidad)
vehiculosb=window(vehiculos,start=1960)
vehiculosbEts=ets(vehiculosb,
model="MAN",
damped=FALSE)
vehiculosbEts2=ets(vehiculosb,
model="MNN",
damped=FALSE)
vehiculosbEts1=accuracy(vehiculosbEts)
vehiculosbEts21=accuracy(vehiculosbEts2)
comparacion=as.data.frame(rbind(vehiculosbEts1,vehiculosbEts21))
rownames(comparacion)=c("Método MAN","Método MNN")
comparacion| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Método MAN | -100.9985 | 908.9398 | 582.1795 | -0.4348065 | 6.089607 | 0.8578568 | 0.1983227 |
| Método MNN | 134.8033 | 928.0030 | 667.6877 | 2.4499440 | 7.437629 | 0.9838555 | 0.2998183 |
- Analizamos el error del modelo MAN (que es con el que nos quedamos) y vemos que no hay ninguna observación importante que se salga de los límites.
error <- residuals(vehiculosbEts)
sderror <- sd(error)
autoplot(error, series="Error",
colour = "black",
xlab = "Semana",
ylab = "Error",
main = "") +
geom_hline(yintercept = c(-3, -2, 2 ,3)*sderror,
colour = c("red", "blue", "blue", "red"), lty = 2) +
scale_x_continuous(breaks= seq(6, 26, 2)) ## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
D) Modelo y las ecuaciones recursivas del mismo
El modelo que tenemos es MAN, es decir:
\(y_{t+h|t}= l_t + hb_t\), con ecuaciones recursivas tal que:
\(l_t=\alpha\,y_t+(1-\alpha)(l_{t-1}+b_{t-1}),\)
\(b_t=\beta(l_t-l_{t-1})+(1-\beta)b_{t-1}.\)
Donde: \(\alpha=0.9998997\),\(\beta=0.1010386\), \(l_t=901.131\) y \(b_t=398.8961\)
El modelo MNN es:
\(y_{t+h|t}=l_t\)
\(l_t=\alpha y_t+(1-\alpha)l_{t-1}\)
Donde: \(\alpha=0.9998995\) y \(l_t=1600.6263\)
tibble(alpha=vehiculosbEts$par["alpha"],beta=vehiculosbEts$par["beta"])| alpha | beta |
|---|---|
| 0.9998997 | 0.1010386 |
tibble(alpha=vehiculosbEts2$par["alpha"])| alpha |
|---|
| 0.9998995 |
E) Validación Cruzada
- El error inicial es 6.69, sin embargo, 5 años más tarde pasa a ser 31.52 aumentando de una forma considerada
k <- 20 #Minimo numero de datos para estimar
h <- 5 #Horizonte de las prediciciones
TT <- length(vehiculos)#Longitud serie
s <- TT - k - h #Total de estimaciones
rmse <- matrix(NA, s + 1, h)
mape <- matrix(NA, s + 1, h)
for (i in 0:s) {
train.set <- subset(vehiculos, start = i + 1, end = i + k)
test.set <- subset(vehiculos, start = i + k + 1, end = i + k + h)
vehiculosbEts=ets(train.set,
model="MAN",
damped=FALSE)
fcast <- forecast(vehiculosbEts, h = h,level=95)
rmse[i + 1,] <- (test.set - fcast$mean)^2
mape[i + 1,] <- 100*abs(test.set - fcast$mean)/test.set
}
rmse.snaive <- sqrt(colMeans(rmse))
mape.snaive <- colMeans(mape)
error=data.frame(t(round(rmse.snaive, 2)))
colnames(error)=c("1A","2A","3A","4A","5A")
mape=data.frame(t(round(mape.snaive, 2)))
colnames(mape)=c("1A","2A","3A","4A","5A")
rownames(mape)=c("Mape")
informacion=rbind(error,mape)
rownames(informacion)=c("Error","Mape")
informacion| 1A | 2A | 3A | 4A | 5A | |
|---|---|---|---|---|---|
| Error | 1054.97 | 1773.37 | 2516.90 | 3545.14 | 4535.10 |
| Mape | 6.69 | 11.67 | 17.71 | 24.62 | 31.52 |
Ejercicio 7: Serie Producción de Chocolate
A) Graficar y comentar las componentes principales
- Esta gráfica es un primer vistazo general de la evolución de la producción de chocolate en Australia
chocolate=read.csv2("Chocolate.csv")
chocolate=ts(chocolate,
start=c(1958,1),
frequency=12)
autoplot(chocolate,
xlab="",
ylab="Chocolate vendido",
main="Chocolate vendido en Australila")- Sin necesidad de crear lineas de marcación, vemos de forma clara como el chocolate ha sido y es un bien cuya producción no ha hecho más que crecer durante el trascurso de toda la serie, en algunos momentos ha sufrido de momentos de decadencia, pero nunca han logrado un efecto dañino para el crecimiento.
chocolateAnual=aggregate(chocolate, FUN=sum)
autoplot(chocolateAnual,
xlab="Años",
ylab="Producción de chocolate",
main="Producción de chocolate en Australia")- Podemos interpretar dicho gráfico como un esquema multiplicativo de la serie, pues a mayor tiempo, mayor valor de los puntos en el gráfico
chocolateDesv=aggregate(chocolate, FUN=sd)
plot(as.numeric(chocolateAnual),as.numeric(chocolateDesv))- En este gráfico de estacionalidad vemos de una forma muy clara como la producción de chocolate aumenta a partir de febrero y comienza a disminuir en septiembre hasta alcanzar el punto más bajo en Enero.
chocolateb=window(chocolate,start=1958)
ggsubseriesplot(chocolateb)+
ylab("Chocolate")- Observando la serie completa vemos que no hay ningún año raro, y que todos parecen bastante corrientes.
ggseasonplot(chocolateb,
year.labels=T,
ylab="Producción de Chocolate")B) Descomposición en regresiones locales
#logchocolate_stl <- stl(log(chocolate),
#s.window = "periodic",
#robust = TRUE)C) Alisado Exponencial
Como hemos visto anteriormente estamos ante un esquema multiplicativo, por lo que tendremos que utilizar un Alisado de Holt-Winters multiplicativo
Este modelo se describe como \(y_{t+h|t}= (l_t +hb_t)*s_{t+h-m(k+1)}\), donde:
\(\alpha=0.1638\) Un valor tan bajo indica que la serie apenas varía con el tiempo
\(\beta=0\) Este valor indica que la pendiente no varía con el paso del tiempo
\(\gamma=0.3358\) Decimos que la estacionalidad permanece más bien constante en el tiempo
chocolateb=window(chocolate,star=1958)
chocolatebEts=ets(chocolateb,
model="MAM",
damped=F)
summary(chocolatebEts)## ETS(M,A,M)
##
## Call:
## ets(y = chocolateb, model = "MAM", damped = F)
##
## Smoothing parameters:
## alpha = 0.1638
## beta = 1e-04
## gamma = 0.3358
##
## Initial states:
## l = 2389.3933
## b = 15.3823
## s = 0.6671 0.8847 0.95 0.9934 1.0893 1.2565
## 1.2352 1.3236 1.1559 1.0285 0.8593 0.5564
##
## sigma: 0.1071
##
## AIC AICc BIC
## 8246.489 8247.926 8316.118
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -15.45195 553.6932 407.4784 -1.031004 8.16777 0.7761905 0.1515466
D) Valores de la última fecha de la serie
- Obtenemos que \(l=9784.753\) y \(b=15.80399\), para poder entender mejor qué es esto, vemos como l es la estimación del nivel de la serie y b es la estimación de la pendiente de la serie.
- Por otro lado \(s_i \forall\,i\in (1,2,...,12)\) tenemos la estimación de la estacionalidad en los meses en orden inverso, s1 es el valor para Diciembre (el último dato), s2 para Noviembre y así hasta s12 que corresponde a Enero
tail(chocolatebEts$states,1)## l b s1 s2 s3 s4 s5
## Dec 1994 9784.753 15.80399 0.9459899 1.056467 1.057233 1.042047 1.140944
## s6 s7 s8 s9 s10 s11 s12
## Dec 1994 0.9891633 0.9695922 0.9650341 0.8533319 1.00468 0.9157171 0.6283261
E) Predicción para los próximos años
chocolatebf=forecast(chocolatebEts,
h=36,
level=95)
chocolatebf$mean## Jan Feb Mar Apr May Jun Jul
## 1995 6157.946 8989.010 9878.183 8403.586 9518.878 9579.161 9788.148
## 1996 6281.023 9168.390 10075.000 8570.762 9707.948 9769.134 9981.965
## 1997 6404.206 9347.926 10271.988 8738.084 9897.181 9959.271 10175.951
## Aug Sep Oct Nov Dec
## 1995 11308.110 10344.393 10511.848 10520.931 9435.683
## 1996 11531.679 10548.593 10719.034 10727.978 9621.088
## 1997 11755.441 10752.970 10926.401 10935.205 9806.655
autoplot(chocolatebf,
xlab="Años",
ylab="Producción de chocolate",
main="Predicción de Chocolate + 3 años",
PI=F)F) Training Set
chocolateIntra <- subset(chocolate, end = length(chocolate) - 120)
chocolateExtra <- subset(chocolate, start = length(chocolate) - 119)
chocolateExtraPre <- snaive(chocolateIntra, h = 56)
accuracy(chocolateExtraPre, chocolateExtra)## ME RMSE MAE MPE MAPE MASE
## Training set 127.3429 538.4386 427.7212 2.443063 10.16616 1.000000
## Test set 1176.8036 1612.5819 1359.9107 16.919180 19.88413 3.179433
## ACF1 Theil's U
## Training set 0.3381947 NA
## Test set 0.5861061 0.6683054