X=read.table(
  "http://freakonometrics.blog.free.fr/public/data/car-sales-quebec.csv",
  header=TRUE,sep=";",nrows=108)
Xt=ts(X[,2],start=c(1960,1),frequency=12)
plot(Xt)

h=12
Y=X[1:(length(Xt)-h),2]
plot(Y,type="l")
Y0=X[(length(Xt)-h+1):(length(Xt)),2]
db=data.frame(Y=Y,X=1:length(Y))
reg=lm(Y~X,data=db)
summary(reg)
abline(reg,col="red")

X=residuals(reg)
plot(X,type="l")

Z=diff(X,12)
plot(Z,type="l")
plot(acf(Z))
plot(pacf(Z))
model1=arima(Z,order=c(12, 0, 0))
model1
plot(acf(residuals(model1)))

model1X=arima(X,order=c(12, 0, 0),
seasonal = list(order = c(0, 1, 0), period = 12))
  
model2X=arima(X,order=c(1, 0, 0),
seasonal = list(order = c(1, 1, 0), period = 12))  
plot(acf(residuals(model2X))) 
Box.test(residuals(model2X),lag=1, type = "Ljung-Box")
Box.test(residuals(model2X),lag=2, type = "Ljung-Box")
Box.test(residuals(model2X),lag=3, type = "Ljung-Box")
Box.test(residuals(model2X),lag=4, type = "Ljung-Box")
Box.test(residuals(model2X),lag=5, type = "Ljung-Box")
Box.test(residuals(model2X),lag=6, type = "Ljung-Box")

plot(acf(X,lag=25))
model3=arima(X,order=c(12, 0, 0))

X1=predict(model1X,12)$pred
X2=predict(model2X,12)$pred
X0=Y0-predict(reg,newdata=data.frame(X=(length(Xt)-h+1):(length(Xt))))
cbind(X0,X1,X2)
plot(X,type="l",xlim=c(0,110))
lines((length(Xt)-h+1):(length(Xt)),X0,type="l",lwd=2)
lines((length(Xt)-h+1):(length(Xt)),X1,type="l",lwd=2,col="red")
lines((length(Xt)-h+1):(length(Xt)),X2,type="l",lwd=2,col="blue")

sum((X0-X1)^2)
sum((X0-X2)^2)

base=read.table(
"http://freakonometrics.blog.free.fr/public/data/car-sales-quebec.csv",
header=TRUE,sep=";",nrows=108)
Y=base[,2]
db=data.frame(Y=Y,X=1:length(Y))
reg=lm(Y~X,data=db)
summary(reg)
plot(Y)
abline(reg,col="red")

X=residuals(reg)
modeltot=arima(X,order=c(1, 0, 0),
seasonal = list(order = c(1, 1, 0), period = 12)) 
Xpred=predict(modeltot,48)
plot(X,type="l",xlim=c(0,150),ylim=c(-10000,10000))
lines(Xpred$pred,col="blue")
lines(Xpred$pred+2*Xpred$se,col="red")
lines(Xpred$pred-2*Xpred$se,col="red")

tendance=predict(reg,newdata=data.frame(X=109:156))
plot(Y,type="l",xlim=c(0,150),ylim=c(5000,30000))
lines(109:156,tendance,col="green")
lines(Xpred$pred+tendance,col="blue")
lines(Xpred$pred+2*Xpred$se+tendance,col="red")
lines(Xpred$pred-2*Xpred$se+tendance,col="red")



report=read.table(
   "C:\\Users\\savsalledecours\\Desktop\\report-headphones.csv",
   skip=4,header=TRUE,sep=",",nrows=464)
 source("http://freakonometrics.blog.free.fr/public/code/H2M.R")
 headphones=H2M(report,lang="FR",type="ts")
 plot(headphones)
 plot(log(headphones))

Y=as.numeric(log(headphones))
plot(Y)
db=data.frame(Y=Y,X=1:length(Y))
reg=lm(Y~X,data=db)
summary(reg)

plot(Y)
plot(acf(Y,25))
plot(pacf(Y,25))
model0=arima(Y,order=c(1, 0, 0))
plot(acf(residuals(model0),37))
model0b=arima(Y,order=c(1, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 12))
model0b
plot(acf(residuals(model0b),37))
Box.test(residuals(model0b),lag=1, type = "Ljung-Box")
Box.test(residuals(model0b),lag=5, type = "Ljung-Box")
Box.test(residuals(model0b),lag=11, type = "Ljung-Box")
Box.test(residuals(model0b),lag=12, type = "Ljung-Box")

model0c=arima(Y,order=c(1, 0, 11),
seasonal = list(order = c(1, 0, 0), period = 12))
model0c
plot(acf(residuals(model0c),37))

abline(reg,col="red")

X=residuals(reg)
plot(X,type="l")
plot(acf(X,37))

Z=diff(X,12)
plot(acf(Z,25))
plot(pacf(Z,25))

model1=arima(Z,order=c(1, 0, 0))
model1
plot(acf(residuals(model1)))

model2=arima(Z,order=c(12, 0, 0))
model2

model2b=arima(Z,order=c(1, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 12))
model2b

plot(acf(X,25))
plot(pacf(X,25))

model3=arima(X,order=c(14, 0, 0))
model3
plot(acf(residuals(model3)))

plot(Y,type="l",xlim=c(90,140),ylim=c(3,4.2))
Y0=predict(model0c,36)$pred
lines(Y0,col="red")
tendance=predict(reg,newdata=data.frame(X=108:143))
lines(108:143,tendance,col="green")

Y3=predict(model3,36)$pred
lines(Y3+tendance,col="blue")
model1c=arima(X,order=c(1, 0, 0),
seasonal = list(order = c(0, 1, 0), period = 12))
model2c=arima(X,order=c(1, 0, 0),
seasonal = list(order = c(1, 1, 0), period = 12))

Y1=predict(model1c,36)$pred
lines(Y1+tendance,col="purple")
Y2=predict(model2c,36)$pred
lines(Y2+tendance,col="orange")

summary(reg)


plot(exp(Y),type="l",xlim=c(90,140),ylim=exp(c(3,4.2)))
plot(exp(Y),type="l",xlim=c(0,140),ylim=c(0,60))
Y0=predict(model0c,36)$pred
lines(exp(.28^2+Y0),col="red")
tendance=predict(reg,newdata=data.frame(X=108:143))
lines(108:143,exp(.28^2+tendance),col="green")

Y3=predict(model3,36)$pred
lines(exp(.28^2+Y3+tendance),col="blue")
model1c=arima(X,order=c(1, 0, 0),
              seasonal = list(order = c(0, 1, 0), period = 12))
model2c=arima(X,order=c(1, 0, 0),
              seasonal = list(order = c(1, 1, 0), period = 12))

Y1=predict(model1c,36)$pred
lines(exp(.28^2+Y1+tendance),col="purple")
Y2=predict(model2c,36)$pred
lines(exp(.28^2+Y2+tendance),col="orange")