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

Machine learning logistic regression for credit modelling in R — R-bloggers

(This article was first published on R Programming – DataScience+, and kindly contributed to R-bloggers) Categories Regression Models Tags Logistic Regression R Programming ROC Machine learning logistic regressions is a widely popular method to model credit modeling. There are excellent and efficient packages in R, that can perform these types of analysis. Typically you will…

via Machine learning logistic regression for credit modelling in R — R-bloggers

5 Alternatives to the Default R Outputs for GLMs and Linear Models — R-bloggers

(This article was first published on R – Displayr, and kindly contributed to R-bloggers) The standard summary outputs from the glm and lm summary methods are a case in point. If you have been using R for as long as I have (19 or 20 years…) you will no doubt have a certain affection for…

via 5 Alternatives to the Default R Outputs for GLMs and Linear Models — R-bloggers

How to add Trend Lines in R Using Plotly — R-bloggers

(This article was first published on R – Displayr, and kindly contributed to R-bloggers) 1. Global trend lines One of the simplest methods to identify trends is to fit a ordinary least squares regression model to the data. The model most people are familiar with is the linear model, but you can add other polynomial…

via How to add Trend Lines in R Using Plotly — R-bloggers

Machine learning at central banks

https://www.bankofengland.co.uk/working-paper/2017/machine-learning-at-central-banks

Published on 01 September 2017

Regression Analysis Essentials For Machine Learning

(This article was first published on Easy Guides, and kindly contributed to R-bloggers) Regression analysis consists of a set of machine learning methods that allow us to predict a continuous outcome variable (y) based on the value of one or multiple predictor variables (x). Briefly, the goal of regression model is to build a mathematical…

via Regression Analysis Essentials For Machine Learning — R-bloggers

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