analisis_cemento.R

In reos156/AnalisisSeriesTiempo:

setwd("D:/GitHub/AnalisisSeriesTiempo")

rm(list=ls())       
graphics.off()

library(forecast)
library(lmtest)
library(xtable)


E = read.table("cementq.dat", header = TRUE)
attach(E)
   
y = ts(E,frequency=4,start=c(1956,1),end=c(1994,3))

#----- validacion cruzada
yi = window(y,start=c(1956,1),end=c(1992,3))
yf = window(y,start=c(1992,4),end=c(1994,3))
m = 8

# alternativa:
# m = 8
# n = length(y)
# yi = ts(y[1:(n-m)],frequency=4)
# yf = ts(y[(n-m+1):n],frequency=4)


lyi = log(yi)
 

#----modelo con tendencia lineal y estacionalidad
#----con variables indicadoras estacionales
require(forecast)

ti = seq(1,length(yi))

It = seasonaldummy(yi)

mod1 = lm(yi ~ ti + It)

summary(mod1)


r1 = mod1$residuals

# analisis de los residuos  

par(mfrow=c(2,3))
plot(ti,r1,type='o',ylab='residuo')
abline(h=0,lty=2)
plot(density(r1),xlab='x',main= '')
qqnorm(r1)
qqline(r1,col=2)
acf(r1,60,ci.type="ma",main="")
pacf(r1,60,main="")


yhat1 = mod1$fitted.values

#----modelo con tendencia lineal y estacionalidad
#----con variables trigonometricas

It.trig = fourier(yi,2)

mod2 = lm(yi ~ ti + It.trig)

summary(mod2)

r2 = mod1$residuals

# analisis de los residuos  

par(mfrow=c(2,2))
plot(ti,r2,type='o',ylab='residuo')
abline(h=0,lty=2)
plot(density(r2),xlab='x',main= '')
qqnorm(r2)
qqline(r2,col=2)

yhat2 = mod2$fitted.values
 
#-----------comparar modelos

source("medidas.r")

M.ind = medidas(mod1,yi,5)
M.trig = medidas(mod2,yi,6)
(M = cbind(M.ind,M.trig))


#---------modelo log cuadratico con estacionalidad
ti2 = ti*ti

mod3 = lm(lyi ~ ti + ti2 + It)

summary(mod3)
 
# El modelo exponencial-cuadratico-estacional
T = length(yi)

Xt = cbind(rep(1,T),ti,ti2,It)

Ds = data.frame(yi,Xt)

theta.0 = mod3$coefficient

mod4 = nls(yi~exp(Xt%*%theta),
data=Ds, start= list(theta=theta.0))

summary(mod4)

yhat4 = fitted(mod4)

M.exp.cuad = medidas(mod4,yi,6)
(M = cbind(M,M.exp.cuad))

#-------- graficas

plot(ti,yi,type='b')
lines(ti,yhat1,col = 'blue')
lines(ti,yhat4,col = 'red')


# -------pronosticos 

T = length(yi)
Itf = seasonaldummy(yi,m)
tf = seq(T+1,T+m,1)
 
pron1 = predict(mod1,
data.frame(t = tf, It=I(Itf)))

# pronosticos con el exp cuad 

tf2 = tf*tf
Xtf = cbind(rep(1,m),tf,tf2,Itf)

pron4 = predict(mod4,
data.frame(Xt = I(Xtf)))

pr = (pron4+pron1)/2

par(mfrow=c(1,1))
plot(tf,yf, type = 'o',ylim=c(1300,2000))
lines(tf,pron1, type = 'b', pch = 3,col='red' )
lines(tf,pron4, type = 'b', pch = 5,col='blue' )

lines(tf,pr, type = 'b', pch = 5,col='orange' )



legend("topleft", 
c("Obs","Indic", "Exp-cuad-estac","HW"), 
pch = c(1, 3, 5,7),
col = c("black","red","blue","brown"))

#-------- graficas
t = seq(1,length(y))

plot(t[100:155],y[100:155],type='b',xlim=c(100,155),lwd=2)
lines(ti[100:147],yhat1[100:147],col = 'blue')
lines(ti[100:147],yhat4[100:147],col = 'red')


lines(tf,pron1, type = 'b', pch = 3,col='red' )
lines(tf,pron4, type = 'b', pch = 5,col='blue' )
lines(tf,pr, type = 'b', pch = 5,col='orange' )

pron4 = ts(pron4,frequency=4,start=c(1992,4),end=c(1994,3))

pr = ts(pr,frequency=4,start=c(1992,4),end=c(1994,3))


A=rbind(
accuracy(pron1,yf), 
accuracy(pron4,yf),
accuracy(pr,yf))

rownames(A) = c("Lin-Ind","Exp-Cuad-Ind","prom")
(A[,c(2,5)])

print(xtable(A))

#------- estimar el efecto de cada trimestre

m.a.1 = lm(yi ~ ti)
Tt = m.a.1$fitted.values
St.et = yi - Tt

It4 = cbind(It,rep(0,T))
It4[,4] = rep(1,T) - apply(It,1,sum)

m.a.2 = lm(St.et ~ It4 - 1)
summary(m.a.2)


#------ filtro brockwell-davis en stl

mod0 = stl(y, s.window = "per" )

mod0 = stl(y, s.window = "per" )

plot(mod0)


#------- Holt Winters
source("medidas.hw.r")

mod5 = HoltWinters(yi, 
seasonal = "additive")

yhat5 = fitted(mod5)[,1]

M.h.w = medidas.hw(mod5,yi,3)
M = cbind(M,M.h.w)
(M)

# pronosticos con holt-winters

pron5 = predict(mod5, m, prediction.interval = TRUE)

lines(tf,pron5[,1], type = 'b', pch = 5,col='brown' )