TRIMESTRALIZACIÓN DEL VAB INDUSTRIAL DE CANTABRIA

Lectura y tratamiento de datos

En indust_m1995.xls están los datos de los indicadores de frecuencia mensual, en indust_T1995.xls los indicadores de frecuencia trimestral, y en CRAC.xls los datos de la serie anual de la contabilidad regional de Cantabria que se quiere trimestralizar.

Datos:

indust_T1995

indust_m1995

CRAC

“`{r}
# Datos mensuales
excel1 <-read.csv(“indust_m1995.csv”, sep=”;”,header = TRUE)
excel3 <-read.csv(“indust_T1995.csv”, sep=”;”,header = TRUE)
excel4<-read.csv(“CRAC.csv”, sep=”;”,header = TRUE)
“`

Creamos objetos ts de las series de Exportaciones (EXPORT), Ocupados (OCUP) en la industria (cifras de la seguridad social), tendencia de la Cartera de Pedidos (TPEDIDOS), Indice de Precios Industriales (IPRI), Indice de producción Industrial (IPI), y Ocupados (EPA) en la Industria en la Encuesta de Población Activa.

“`{r}
EXPORT=ts(na.omit(excel1$EXPORT),frequency = 12,start=c(1995,1))
OCUP=ts(na.omit(excel1$OCUP),frequency = 12,start=c(1995,1))
TPEDIDO=ts(na.omit(excel1$TPEDIDO),frequency = 12,start=c(1995,1))
IPRI=ts(na.omit(excel1$IPRI),frequency = 12,start=c(1995,1))
IPI=ts(na.omit(excel1$IPI),frequency = 12,start=c(1995,1))
EPA=ts(na.omit(excel3$EPA),frequency = 4,start=c(1995,1))
“`

Convertimos las series de indicadores de frecuencia mensual a trimestral:

“`{r}
EXPORT=aggregate.ts(EXPORT,nfrequency = 4,mean)
OCUP=aggregate.ts(OCUP,nfrequency = 4,mean)
TPEDIDO=aggregate.ts(TPEDIDO,nfrequency = 4,mean)
IPRI=aggregate.ts(IPRI,nfrequency = 4,mean)
IPI=aggregate.ts(IPI,nfrequency = 4,mean)
“`

Deflactamos la serie de exportaciones por el índice de precios industriales.

“`{r}
EXPORT=EXPORT*100/IPRI
“`

Desestacionalizamos

Obtenemos las series de tendencia, estacionalidad con descomponer, con modelo aditivo.

“`{r}
library(descomponer)
#EXPORT
EXPORT.s=ts(descomponer(EXPORT,4,1)$datos$TDST,frequency=4,start=c(1995,1))
EXPORT.sa=ts(descomponer(EXPORT,4,1)$datos$ST,frequency=4,start=c(1995,1))
EXPORT.t=ts(descomponer(EXPORT,4,1)$datos$TD,frequency=4,start=c(1995,1))
plot(EXPORT)
lines(EXPORT.s,col=”red”)
lines(EXPORT.t,col=”blue”)
#OCUP
OCUP.s=ts(descomponer(OCUP,4,1)$datos$TDST,frequency=4,start=c(1995,1))
OCUP.sa=ts(descomponer(OCUP,4,1)$datos$ST,frequency=4,start=c(1995,1))
OCUP.t=ts(descomponer(OCUP,4,1)$datos$TD,frequency=4,start=c(1995,1))
plot(OCUP)
lines(OCUP.s,col=”red”)
lines(OCUP.t,col=”blue”)
#TPEDIDO
TPEDIDO.s=ts(descomponer(TPEDIDO,4,1)$datos$TDST,frequency=4,start=c(1995,1))
TPEDIDO.sa=ts(descomponer(TPEDIDO,4,1)$datos$ST,frequency=4,start=c(1995,1))
TPEDIDO.t=ts(descomponer(TPEDIDO,4,1)$datos$TD,frequency=4,start=c(1995,1))
plot(TPEDIDO)
lines(TPEDIDO.s,col=”red”)
lines(TPEDIDO.t,col=”blue”)
#IPI
IPI.s=ts(descomponer(IPI,4,1)$datos$TDST,frequency=4,start=c(1995,1))
IPI.sa=ts(descomponer(IPI,4,1)$datos$ST,frequency=4,start=c(1995,1))
IPI.t=ts(descomponer(IPI,4,1)$datos$TD,frequency=4,start=c(1995,1))
plot(IPI)
lines(IPI.s,col=”red”)
lines(IPI.t,col=”blue”)
#EPA
EPA.s=ts(descomponer(EPA,4,1)$datos$TDST,frequency=4,start=c(1995,1))
EPA.sa=ts(descomponer(EPA,4,1)$datos$ST,frequency=4,start=c(1995,1))
EPA.t=ts(descomponer(EPA,4,1)$datos$TD,frequency=4,start=c(1995,1))
plot(EPA)
lines(EPA.s,col=”red”)
lines(EPA.t,col=”blue”)
“`

Pronósticos de la serie desestacionalizada

Con tslm de forecast pronosticamos: tendencias y estacionalidades. Prolongamos la serie original (o) con el pronostico de tendencia y estacionalidad, la serie de tendencia y estacionalidad (f) con el pronostico de tendencia y estacionalidad, y la serie de tendencia (t) con el pronóstico de tendencia.

“`{r}
library(forecast)
#EXPORT
fit1 <- auto.arima(EXPORT.sa)
fit2 <- auto.arima(EXPORT.t)
plot(forecast(fit1, h=20))
plot(forecast(fit2, h=20))
fit3 = forecast(fit1, h=20)
fit4 = forecast(fit2, h=20)
EXPORT.f=ts(c(EXPORT.s,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
EXPORT.o=ts(c(EXPORT,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
EXPORT.t=ts(c(EXPORT.t,fit4$mean),frequency = 4,start=c(1995,1))
plot(EXPORT.o)
lines(EXPORT.f,col=”blue”)
lines(EXPORT.t,col=”red”)
#OCUP
fit1 <- auto.arima(OCUP.sa)
fit2 <- auto.arima(OCUP.t)
plot(forecast(fit1, h=20))
plot(forecast(fit2, h=20))
fit3 = forecast(fit1, h=20)
fit4 = forecast(fit2, h=20)
OCUP.f=ts(c(OCUP.s,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
OCUP.o=ts(c(OCUP,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
OCUP.t=ts(c(OCUP.t,fit4$mean),frequency = 4,start=c(1995,1))
plot(OCUP.o)
lines(OCUP.f,col=”blue”)
lines(OCUP.t,col=”red”)
# TPEDIDO
fit1 <- auto.arima(TPEDIDO.sa)
fit2 <- auto.arima(TPEDIDO.t)
plot(forecast(fit1, h=20))
fit3 = forecast(fit1, h=20)
fit4 = forecast(fit2, h=20)
TPEDIDO.f=ts(c(TPEDIDO.s,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
TPEDIDO.o=ts(c(TPEDIDO,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
TPEDIDO.t=ts(c(TPEDIDO.t,fit4$mean),frequency = 4,start=c(1995,1))
plot(TPEDIDO)
lines(TPEDIDO.f,col=”blue”)
lines(TPEDIDO.t,col=”red”)
#IPI
fit1 <- auto.arima(IPI.sa)
fit2 <- auto.arima(IPI.t)
plot(forecast(fit1, h=20))
fit3 = forecast(fit1, h=20)
fit4 = forecast(fit2, h=20)
IPI.f=ts(c(IPI.s,c(fit3$mean)+c(fit4$mean)),frequency = 4,start=c(1995,1))
IPI.o=ts(c(IPI,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
IPI.t=ts(c(IPI.t,fit4$mean),frequency = 4,start=c(1995,1))
plot(IPI)
lines(IPI.f,col=”blue”)
lines(IPI.t,col=”red”)
# EPA
fit1 <- auto.arima(EPA.sa)
fit2 <- auto.arima(EPA.t)
plot(forecast(fit1, h=20))
fit3 = forecast(fit1, h=20)
fit4 = forecast(fit2, h=20)
EPA.f=ts(c(EPA.s,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
EPA.o=ts(c(EPA,fit3$mean+fit4$mean),frequency = 4,start=c(1995,1))
EPA.t=ts(c(EPA.t,fit4$mean),frequency = 4,start=c(1995,1))
plot(EPA)
lines(EPA.f,col=”blue”)
lines(EPA.t,col=”red”)

“`

Agregación de series

Obtenemos los indicadores en frecuencias anuales desde las frecuencias trimestrales.

“`{r}
EXPORT.A=aggregate.ts(EXPORT.f,nfrequency = 1,mean)
OCUP.A=aggregate.ts(OCUP.f,nfrequency = 1,mean)
TPEDIDO.A=aggregate.ts(TPEDIDO.f,nfrequency = 1,mean)
IPI.A=aggregate.ts(IPI.f,nfrequency = 1,mean)
EPA.A=aggregate.ts(EPA.f,nfrequency = 1,mean)
“`

Creamos varios data frame con las series anuales para el periodo 2000-2019. En cada data frame ha

“`{r}
fin=2019
# SELECCION BASE
INDUS_CAN.A=data.frame(EXPORT=window(EXPORT.A,start=c(2000),end=c(fin)),OCUP=window(OCUP.A,start=c(2000),end=c(fin)),IPI=window(IPI.A,start=c(2000),end=c(fin)))
INDUS_CAN.A.ts=ts(INDUS_CAN.A,start=c(2000))
INDUS_CAN.A.tasa=diff(INDUS_CAN.A.ts)/INDUS_CAN.A.ts
colnames(INDUS_CAN.A.tasa)=c(“EXPORT”,”OCUP”,”IPI”)
# SELECCION 1
INDUS_CAN.A.1=data.frame(EXPORT=window(EXPORT.A,start=c(2000),end=c(fin)),OCUP=window(OCUP.A,start=c(2000),end=c(fin)),TPEDIDO=window(TPEDIDO.A,start=c(2000),end=c(fin)),IPI=window(IPI.A,start=c(2000),end=c(fin)))
INDUS_CAN.A.1.ts=ts(INDUS_CAN.A.1,start=c(2000))
INDUS_CAN.A.1.tasa=diff(INDUS_CAN.A.1.ts)/INDUS_CAN.A.1.ts
colnames(INDUS_CAN.A.1.tasa)=c(“EXPORT”,”OCUP”,”TPEDIDO”,”IPI”)
# SELECCION 2
INDUS_CAN.A.2=data.frame(EXPORT=window(EXPORT.A,start=c(2000),end=c(fin)),IPI=window(IPI.A,start=c(2000),end=c(fin)),EPA=window(EPA.A,start=c(2000),end=c(fin)))
INDUS_CAN.A.2.ts=ts(INDUS_CAN.A.2,start=c(2000))
INDUS_CAN.A.2.tasa=diff(INDUS_CAN.A.2.ts)/INDUS_CAN.A.2.ts
colnames(INDUS_CAN.A.2.tasa)=c(“EXPORT”,”IPI”,”EPA”)
# SELECCION 3
INDUS_CAN.A.3=data.frame(OCUP=window(OCUP.A,start=c(2000),end=c(fin)),IPI=window(IPI.A,start=c(2000),end=c(fin)))
INDUS_CAN.A.3.ts=ts(INDUS_CAN.A.3,start=c(2000))
INDUS_CAN.A.3.tasa=diff(INDUS_CAN.A.3.ts)/INDUS_CAN.A.3.ts
colnames(INDUS_CAN.A.3.tasa)=c(“OCUP”,”IPI”)
“`

Leemos la serie anual del Indice de Volumen del VAB de la industria de la CRE y obtenemos la tasas de crecimento anual:

“`{r}
VAB=ts(na.omit(excel4$VAB2_Co),start=c(2000))
VAB.tasa=diff(VAB)/VAB
“`

Índices con Análisis factorial

Elaboramos un indicador de síntesis con la selección básica de indicadores, utilizando la librería “rela”, se correlaciona el indicador obtenido con la tasa anual del crecimiento del VAB en la industria, y se proyecta la serie anual del VAB hasta 2019.

“`{r}
library(rela)
res <- paf(as.matrix(INDUS_CAN.A.tasa))
summary(res) # Obtenemos el KMO con MSA, para determinar el número de factores
# calculamos el chi-square del Bartlett’s sphericity test, la comunalidad y
# pesos de los factores.
barplot(res$Eigenvalues[,1]) # Representamos los eigenvalues.
resv <- varimax(res$Factor.Loadings) # Rotación Varimax
print(resv)
#barplot(sort(colSums(loadings(resv)^2),decreasing=TRUE)) # Gráfico utilizando los pesos rotados.
scores <- scale(as.matrix(INDUS_CAN.A.tasa)) %*% as.matrix(resv) # Pesos factoriales, para obtener los indices tambien se puede
# Intervenciones
d13=c(rep(0,12),1,rep(0,(17-13)))
d14=c(rep(0,13),1,rep(0,(17-14)))
# Modelo 1: Estimacion MCO
indice1=ts(scores[,1],start=c(2001))
indice1.1=window(indice1,end=c(2017))
mod1=lm(VAB.tasa~0+indice1.1)
summary(mod1)
plot(mod1)
# Grafico
plot(VAB.tasa)
lines(ts(mod1$fitted,start=c(2001)),col=”red”)
# Pronósticos
indice1.1=window(indice1,start=c(2018),end=c(2019))
nuevos=data.frame(indice1.1)
VAB.tasa.fa=ts(predict(mod1,nuevos,interval=”prediction”),start=c(2018),end=c(2019))
VAB.tasa.fa
VAB.2018=window(VAB,start=c(2017),end=c(2017))*c((1+window(VAB.tasa.fa[,”fit”],start=c(2018),end=c(2018))))
VAB.2019=VAB.2018*c((1+window(VAB.tasa.fa[,”fit”],start=c(2019),end=c(2019))))
VAB1=ts(c(VAB,VAB.2018,VAB.2019),start=c(2000))
plot(VAB1)
“`

Índices con el método de Granger y Newbold

Elaboramos un indicador de síntesis para la industria utilizando el método de Granger y Newbold y la selección básica de indicadores, se correlaciona el indicador obtenido con la tasa anual del crecimiento del VAB en la industria, y se proyecta la serie anual del VAB hasta 2019 :

“`{r}

# Indice
sd1=1/sd(lm(VAB.tasa~window(INDUS_CAN.A.tasa[,”EXPORT”],end=c(2017)))$residuals)
sd2=1/sd(lm(VAB.tasa~window(INDUS_CAN.A.tasa[,”OCUP”],end=c(2017)))$residuals)
sd3=1/sd(lm(VAB.tasa~window(INDUS_CAN.A.tasa[,”IPI”],end=c(2017)))$residuals)
sdT=sd1+sd2+sd3
indice2=(sd1/sdT)*INDUS_CAN.A.tasa[,”EXPORT”]+(sd2/sdT)*INDUS_CAN.A.tasa[,”OCUP”]+(sd3/sdT)*INDUS_CAN.A.tasa[,”IPI”]
indice2=ts(indice2,start=c(2001))
indice2.1=window(indice2,end=c(2017))
# Estimacion MCO
mod2=lm(VAB.tasa~0+indice2.1)
#mod2=lm(VAB.tasa~0+indice2.1)
summary(mod2)
plot(mod2)
# Grafico
plot(VAB.tasa)
lines(ts(mod2$fitted,start=c(2001)),col=”red”)
# Pronósticos
indice2.1=window(indice2,start=c(2018),end=c(2019))
nuevos=data.frame(indice2.1)
VAB.tasa.f=ts(predict(mod2,nuevos,interval=”prediction”),start=c(2018),end=c(2019))
VAB.tasa.f
VAB.2018=window(VAB,start=c(2017),end=c(2017))*c((1+window(VAB.tasa.f[,”fit”],start=c(2018),end=c(2018))))
VAB.2019=VAB.2018*c((1+window(VAB.tasa.f[,”fit”],start=c(2019),end=c(2019))))
VAB2=ts(c(VAB,VAB.2018,VAB.2019),start=c(2000))
plot(VAB2)
“`

Trimestralizar

Trimestralizamos la serie original con el método de Granger y Newbold, utilizando la librería tempdisagg, y el método de Chow-Lee.

“`{r}

library(tempdisagg)
indice.o=(sd1/sdT)*scale(window(EXPORT.o,start=c(2000,1),end=c(2018,1)))+(sd2/sdT)*scale(window(OCUP.o,start=c(2000,1),end=c(2018,1)))+(sd3/sdT)*scale(window(IPI.o,start=c(2000,1),end=c(2018,1)))
indice.o=ts(indice.o,frequency=4,start=c(2000,1))
mod.o=td(VAB1~indice.o,method = “chow-lin-maxlog”)
summary(mod.o)
plot(mod.o)
resultado.o=predict(mod.o)
resultado.o
plot(resultado.o)
“`

Trimestralizamos la serie de tendencia, y la serie de tendencia y estacionalidad, utilizando la libreria tempdisagg, y el método de Chow-Lee, representamos los resultados finales:

“`{r}
indice.f=(sd1/sdT)*scale(window(EXPORT.f,start=c(2000,1),end=c(2018,1)))+(sd2/sdT)*scale(window(OCUP.f,start=c(2000,1),end=c(2018,1)))+(sd3/sdT)*scale(window(IPI.f,start=c(2000,1),end=c(2018,1)))
indice.f=ts(indice.f,frequency=4,start=c(2000,1))
mod.f=td(VAB1~indice.f,method = “chow-lin-maxlog”)
summary(mod.f)
plot(mod.f)
indice.T=(1/sd1)/(1/sdT)*scale(window(EXPORT.t,start=c(2000,1),end=c(2018,1)))+(1/sd2)/(1/sdT)*scale(window(OCUP.t,start=c(2000,1),end=c(2018,1)))+(1/sd3)/(1/sdT)*scale(window(IPI.t,start=c(2000,1),end=c(2018,1)))
indice.T=ts(indice.T,frequency=4,start=c(2000,1))
mod.T=td(VAB1~indice.T,method = “chow-lin-maxlog”)
summary(mod.T)
plot(mod.T)
resultado.f=predict(mod.f)
resultado.f
resultado.T=predict(mod.T)
resultado.T
plot(resultado.o)
lines(resultado.T,col=”red”)
lines(resultado.f,col=”blue”)
“`

Resutados en word

Trimestralizacion industria

En R pub:

https://rpubs.com/PacoParra/443765

 

 

 

Anuncios

Modelizar el consumo eléctrico

A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à ca plot(electricite[passe,”Load”],type=”l”) On peut tenter un modèle assez simple, où la consommation à la date Y_t est fonction d’une tendance linéaire a+bt, de la position dans l’année (sous une forme…

via Modéliser la consommation électrique — Freakonometrics

ESS Guidelines on Seasonal Adjustment. Edición 2015

Eurostat ha publicat l’edició 2015 de les recomanacions per desestacionalització: ESS guidelines on seasonal adjustment “The revised ESS Guidelines on Seasonal Adjustment present both theoretical aspects and practical implementation issues in a friendly and easy to read framework, thereby addressing both experts and non-experts in seasonal adjustment. They meet the requirement of principle 7 (Sound […]

via ESS Guidelines on Seasonal Adjustment. Edició 2015 — Bloc d’estadística oficial

Descomposición temporal con R

 

Resumen:

Evaluación de diferentes funciones descomposición temporal de R. Se utiliza como ejemplo la serie mensual de las concentraciones atmosféricas de CO2 en partes por millón (ppm) en Mauna Loa (Hawai),desde 1959 a 1979 .

Las funciones de R analizadas son: decompose, stl, decomp y descomponer. En los análisis gráficos muestran una gran simulitud, tanto en lo relativo a la serie de tendencia T(t), como a la estacionalidad S(t). La serie irregular I(t) acepta la hipótesis de normalidad en los test estadísticos KS y CVM en todos los casos, pero con autocorrelación serial.

 

descomposición temporal con R

 

En R-Pub:

Descomposición temporal con R

timekit: Time Series Forecast Applications Using Data Mining

(This article was first published on business-science.io – Articles, and kindly contributed to R-bloggers) The timekit package contains a collection of tools for working with time series in R. There’s a number of benefits. One of the biggest is the ability to use a time series signature to predict future values (forecast) through data mining…

via timekit: Time Series Forecast Applications Using Data Mining — R-bloggers