Series de tiempo en R

Series de tiempo: análisis de tendencia, descomposición y modelos ARIMA y redes neuronales,


###consumo aparente de vehículos(2000-2015) fuente:ANDI
ti<-seq(2000,2015,1);ti
consumo<-c(61834,67525,98702,75605,104623,150421,225139,273367,213668,177976,267472,351012,325278,
294857,339670,280539)
plot(ti,consumo,type="o")

##tendencia
##modelo lineal  y=a+bt
t<-c(seq(1,16,1))
plot(t,consumo,type="o")
ml<-lm(consumo~t);ml
summary(ml)
pre1<-predict(ml);pre1
lines(pre1,lwd=2,col="red")
res1<-resid(ml);res1
cbind(consumo,pre1,res1)
nueva1 <- data.frame(t = seq(17,20))
predict(ml,nueva1)
AIC(ml)

nuevo1 = seq(17,20,1)
pr1 = predict(ml,data.frame(t=nuevo1));pr1


##modelo cuadrático  y=a+bt+bt^2
t2<-t^2
mc<-lm(consumo~t+t2);mc
summary(mc)
plot(t,consumo,type="o")
pre2<-predict(mc);pre2
lines(pre2,lwd=2,col="red")
res2<-resid(mc);res2
cbind(consumo,pre2,res2)
AIC(mc)
nueva2 <- data.frame(t = seq(17,20),t2)
predict(mc,nueva2)

##modelo cubico
t3<-t^3
mcu<-lm(consumo~t+t2+t3);mcu
summary(mcu)
plot(t,consumo,type="o")
pre3<-predict(mcu);pre3
lines(pre3,lwd=2,col="red")
res3<-resid(mcu);res3
cbind(consumo,pre3,res3)
AIC(mcu)
nueva3 <- data.frame(t = seq(17,20),t3)
predict(ml,nueva1)

##exponencial y=a*exp(bt)
m4<-nls(consumo~a*exp(b*t),start=list(a=20000,b=0.1));m4
summary(m4)
plot(t,consumo,type="o")
lines(pre4,lwd=2,col="red")
pre4<-predict(m4);pre4
res4<-resid(m4)
cbind(consumo,pre4,res4)
AIC(m4)

nuevo4 = seq(17,20,1)
pr4 = predict(m4,data.frame(t=nuevo4));pr4


##Reciproco y=a+b/t
m5<-nls(consumo~a+b*1/t,start=list(a=10000,b=0.1));m5
summary(m5)
plot(t,consumo,type="o")
pre5<-predict(m5);pre5
lines(pre5,lwd=2,col="red")
res5<-resid(m5)
cbind(consumo,pre5,res5)
AIC(m5)

nuevo5 = seq(17,20,1)
pr5 = predict(m5,data.frame(t=nuevo5));pr5

### exponencial y=a*t^b
m6<-nls(consumo~a*t^b,start=list(a=10000,b=0.1));m6
summary(m6)
plot(t,consumo,type="o")
pre6<-predict(m6);pre6
lines(pre6,lwd=2,col="red")
res6<-resid(m6)
cbind(consumo,pre6,res6)
AIC(m6)

nuevo6 = seq(17,20,1)
pr6 = predict(m6,data.frame(t=nuevo6));pr6


##exponencia y=exp(a+bt)
m7<-nls(consumo~exp(a+b*t),start=list(a=20,b=1));m7
summary(m7)
plot(t,consumo,type="o")
pre7<-predict(m7);pre7
lines(pre7,lwd=2,col="red")
res7<-resid(m7)
cbind(consumo,pre7,res7)
AIC(m7)

nuevo7 = seq(17,20,1)
pr7 = predict(m7,data.frame(t=nuevo7));pr7

##

###
library(nnet)
library(forecast)
library(caret)
library(neuralnet)

plot(consumo,type="l")
aj1=auto.arima(consumo)
aj1
summary(aj1)
pre20<-predict(aj1);pre20
plot(forecast(aj1,h=5))
points(1:length(consumo),fitted(aj1),type="l",col="green")

aj2 <- ets(consumo)
summary(aj2)
pre21<-predict(aj2);pre21
plot(forecast(aj2,h=5))
points(1:length(consumo),fitted(aj2),type="l",col="red")

aj3 <- nnetar(consumo)
summary(aj3)
predict(aj3)
plot(forecast(aj3,h=5))
points(1:length(consumo),fitted(aj3),type="l",col="green")


#####
####

2010 2011 2012 2013 2014
##sacrificio de ganado
library(TTR)

datos<-c(295839,189120,209502,339964,354888,279324,190823,215406,307956,314295,
297428,215370,241733,307595,326528,288924,196892,218068,350029,309083,
297811,226628,252026,362010,347737,297639,235449,253180,338591,317963,
314819,229895,247540,346190,349251,298125,236100,256652,341287,331098,
299129,232521,237528,333125,326818,311325,233826,261270,358702,339998,
307341,237653,261780,344347,315589,335957,294523,321571,356546,349263)
y = ts(datos,frequency=12,start=c(2010,01))
y
plot.ts(y)
suav<-SMA(y,n=12)
lines(suav,col="red",lwd=2)
des<-decompose(y)
des
plot(des)
Tend<-decompose(y)$trend
plot(Tend)

ad<-y-des$seasonal
plot.ts(y)
lines(ad,col="red",lwd=2)
plot(ad)
HW<-HoltWinters(y, beta=FALSE, gamma=FALSE)
HW
HW1<-HoltWinters(y, gamma=FALSE)
HW1
HW2<-HoltWinters(y);HW2
fitted(HW1)
plot(HW1)

fore1<-forecast.HoltWinters(HW1, h=12);fore1
plot.forecast(fore1)
acf(fore1$residuals, lag.max=20)
Box.test(fore1$residuals, lag=20, type="Ljung-Box")
plot.ts(fore1$residuals)

di <- diff(y, differences=1);di
plot.ts(di)

di2<- diff(y, differences=2);di2
plot.ts(di2)
acf(di2, lag.max=10)
pacf(di2, lag.max=10)

ari<-arima(y, order=c(0,1,1));ari
forecast(ari)
tsdiag(ari)
accuracy(ari)

ajus <- ets(y);ajus
forecast(ajus,1)
plot(forecast(ajus, 5))

seasonplot(y)

plot(y,type="l")
aj5=auto.arima(y)
aj5
summary(aj5)
pre50<-forecast(aj5);pre50
plot(forecast(aj5,h=12))
points(1:length(y),fitted(aj5),type="l",col="green")

aj6 <- ets(y)
summary(aj6)
pre51<-predict(aj6);pre51
plot(forecast(aj6,h=12))
points(1:length(y),fitted(aj6),type="l",col="red")

aj7 <- nnetar(y)
summary(aj7)
predict(aj7)
plot(forecast(aj7,h=12))
points(1:length(y),fitted(aj7),type="l",col="green")

##modelo AR
mo1<-ar(y);mo1
pro1<-forecast(mo1);pro1
plot(y,type="l")
lines(mo1)

plot(pro1)
pro2<-predict(mo1,n.ahead=12);pro2

mo2<-ma(y,order=2);mo2
plot(y,type="l")
lines(mo2,col="red")
pred<-predict(mo2);pred
mo3<-arima(y,order=c(2,2,0))
mo3
predict(mo3)