#IGE, Ministerio de Fomento. Estatística de construción de edificios
#Datos definitivos dende o ano 1990 ata o ano 2008


setwd("C:/Users/leyenda/Desktop/master/Series_de_tiempo")

edificios <- scan("edificios.txt") #importa las observaciones
edificios.ts <- ts(edificios, start=c(1990, 1), frequency=12)


# Reservamos los datos del último año para valorar el comportamiento del modelo Box-Jenkins que ajustaremos al resto de datos.

edificios1.ts <- window(edificios.ts, end=c(2007,12))
edificios2.ts <- window(edificios.ts, start=c(2008,1))

# Gráfico de secuencia
windows()
plot(edificios1.ts, type="o", xlab="Año", ylab="Edifios construídos")

# Fas y fap estimadas
windows()
par(mfrow=c(2,1))
acf(edificios1.ts, lag.max=length(edificios1.ts)/4, main="")
pacf(edificios1.ts, lag.max=length(edificios1.ts)/4, main="")

# Presencia de tendencia: diferenciación regular
edificios1.ts.dif <- diff(edificios1.ts,1)

windows()
plot(edificios1.ts.dif, type="o", xlab="Mes", ylab="")

windows()
par(mfrow=c(2,1))
acf(edificios1.ts.dif, xlab="Retardo", main="", lag.max=length(edificios1.ts)/4)
pacf(edificios1.ts.dif, xlab="Retardo", main="", lag.max=length(edificios1.ts)/4)



### CUIDADO !!!: la función mejor.arma localiza el ARMA que minimiza uno de los 3 criterios (AIC, AICC o BIC); NO localiza ARIMAs.
###          Por tanto, si lo que buscamos es un ARIMA (d o D no nulo/s), debemos introducir como x la serie convenientemente diferenciada.
 
mejor.arma(x=edificios1.ts.dif, ord.max=c(3,3,2,2), incluir.media="SI", metodo=NULL, criterio="BIC", dist.max.crit=2)

##### Criterio BIC; Método de estimación: mínimos cuadrados condicionados

mejor.arma(x=log.electricidad1.ts.dif.dif12, ord.max=c(3,3,2,2), incluir.media="SI", metodo="CSS", criterio="BIC", dist.max.crit=2)

### Ajuste por máxima verosimilitud

ajuste <- arima(log.electricidad1.ts.dif.dif12, order=c(1,0,2), seasonal=list(order=c(0,0,1), period=12), include.mean=TRUE, optim.control=list(maxit=500))
ajuste

# Media/cte no significativa. Volvemos al proceso de selección, ahora entre modelos sin media/cte (media/cte=0)

mejor.arma(x=log.electricidad1.ts.dif.dif12, ord.max=c(3,3,2,2), incluir.media="NO", metodo="CSS", criterio="BIC", dist.max.crit=2)


###### Modelos tentativos (BIC): ARMA(1,2) x ARMA(0,1)_12 y ARMA(1,3) x ARMA(0,1)_12, ambos sin media/cte y para los (log) datos diferenciados; 
###### esto es, ARIMA(1,1,2) x ARMA(0,1,1)_12 y ARIMA(1,1,3) x ARMA(0,1,1)_12 sin cte para los (log) datos sin diferenciar.

### Ajustes por máxima verosimilitud
ajuste1 <- arima(log.electricidad1.ts, order=c(1,1,2), seasonal=list(order=c(0,1,1), period=12), optim.control=list(maxit=500))
ajuste1

ajuste2 <- arima(log.electricidad1.ts, order=c(1,1,3), seasonal=list(order=c(0,1,1), period=12), optim.control=list(maxit=500))
ajuste2

AIC(ajuste1, ajuste2, k=log(length(log.electricidad1.ts.dif.dif12)))

### El primero tiene un valor BIC (MV) claramente inferior (y menos parámetros). Si pasa el análisis de residuos, nos quedaremos con él

###### Modelo tentativo (BIC): ARMA(1,2) x ARMA(0,1)_12 sin media/cte para los (log) datos diferenciados; 
###### esto es, ARIMA(1,1,2) x ARMA(0,1,1)_12 para los (log) datos sin diferenciar.

### Ajuste por máxima verosimilitud

ajuste <- ajuste1
ajuste

# Todos los coeficientes son significativos 

### Análisis de residuos

# Gráfico de secuencia de los residuos del ajuste arma anterior
windows()
plot(ajuste$resid, type="o", xlab="Tiempo", ylab="Residuos")
abline(0,0)

# Gráfico Q-Q normal
windows()
qqnorm(ajuste$resid)

# Fas estimadas y Ljung-Box
windows()
tsdiag(ajuste, gof.lag=25)

# mu_a = 0
t.test(ajuste$resid, mu=0)

# Normalidad
library(tseries)
jarque.bera.test(ajuste$resid)
shapiro.test(ajuste$resid)

## Conclusión: El modelo ajustado ARIMA(1,1,2) x ARMA(0,1,1)_12 (sin cte (0)
##             puede ser utilizado como generador de la serie del (log) Consumo de Electricidad. 
##             Las innovaciones se pueden considerar son gaussianas.


############################## Predicción
############ Utilizando el modelo ajustado ARIMA(1,1,2) x ARMA(0,1,1)_12, predecimos el consumo del último año (que no fue utilizado para construir y ajustar el modelo)

log.electricidad1.pr <- predict(ajuste, n.ahead=12)

# Predicciones puntuales (del log-consumo)
log.electricidad1.pr$pred

# Desviación típica de las predicciones (del log-consumo)
log.electricidad1.pr$se

# Intervalos de predicción (del log-consumo)
# Extremo inferior
log.electricidad1.pr$pred - 1.96*log.electricidad1.pr$se

# Extremo superior
log.electricidad1.pr$pred + 1.96*log.electricidad1.pr$se

# Intervalos de predicción (del consumo)
ext.inf <- exp(log.electricidad1.pr$pred - 1.96*log.electricidad1.pr$se)
ext.sup <- exp(log.electricidad1.pr$pred + 1.96*log.electricidad1.pr$se)

ext.inf
ext.sup

ts.plot(cbind(window(electricidad1.ts, start=2000), electricidad2.ts, exp(log.electricidad1.pr$pred), ext.inf, ext.sup), col=c("black", "green", "blue", "red", "red"), lty=c(1,1,1,2,2), type="o", xlab="Mes", ylab="Consumo de electricidad")
legend(x=2000, y=max(ext.sup), legend = c("Serie temporal", "Valores a predecir", "Predicciones", "Intervalos de predicción"), col=c("black", "green", "blue", "red"), lty=c(1,1,1,2), pch=1, bty="n")