###################################################################
#                         Ejemplo Shumway                         #
###################################################################

###############################
#Se establece el directorio de
#trabajo y se carga la información 
#de muertes cardiovasculares
###############################


setwd("C:/Users/Mayela/Desktop/Expo_series")
library(astsa)
require(astsa)
require(lmtest)
data(cmort)
data(tempr)
data(part)

par(mfrow=c(3,1))
plot(cmort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(tempr, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
dev.new() # open a new graphic device for the scatterplot matrix
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part))
temp = tempr-mean(tempr) # center temperature
temp2 = temp^2
trend = time(cmort) # time
#regresion
fit = lm(cmort~ trend + temp + temp2 + part, na.action=NULL)

#plot, acf y pacf de los residuales
#se observa que los residuales son un AR(2)
res<-fit$residuals
plot(res,main="Residuales")
acf(res,main="ACF de los residuales")
pacf(res,main="PACF de los residuales")

#prueba de Durbin Watspon
dwtest(cmort~ trend + temp + temp2 + part)


#Se ajusta el modelo de regresión a los datos 
(fit2 <- arima(cmort, order=c(2,0,0),xreg=cbind(trend,temp,temp2
                                              ,part)))

#plot, acf y pacf de los residuales finales
#despues de cconsiderar que los errores son AR(2)
plot(fit2$residuals, ylab="res 2",main="Residuales finales")
acf(fit2$residuals,main="ACF de Residuales finales")
pacf(fit2$residuals ,main="PACF de Residuales finales")

#prueba de Durbin Watson
coef<-fit2$coef
n<-length(cmort)
cmortnew<-cmort[-c(1,2)]-coef[1]*cmort[-c(1,n)]-coef[2]*cmort[-c(n,n-1)]
tempnew<-temp[-c(1,2)]-coef[1]*temp[-c(1,n)]-coef[2]*temp[-c(n,n-1)]
temp2new<-temp2[-c(1,2)]-coef[1]*temp2[-c(1,n)]-coef[2]*temp2[-c(n,n-1)]
trendnew<-trend[-c(1,2)]-coef[1]*trend[-c(1,n)]-coef[2]*trend[-c(n,n-1)]
partnew<-part[-c(1,2)]-coef[1]*part[-c(1,n)]-coef[2]*part[-c(n,n-1)]
fitnew = lm(cmortnew~ trendnew + tempnew + temp2new + partnew, na.action=NULL)
dwtest(cmortnew~ trendnew + tempnew + temp2new + partnew)










#regresión para estimar los parametros del AR(2)

#fitres<-lm(res[-c(1,2)]~0+res[-c(1,n)]+res[-c(n,n-1)])
#plot(fitres$residuals)
#coef<-fitres$coefficients


#transformar el modelo para que tenga residuales indepndientes
#cmortnew<-cmort[-c(1,2)]-coef[1]*cmort[-c(1,n)]-coef[2]*cmort[-c(n,n-1)]
#tempnew<-temp[-c(1,2)]-coef[1]*temp[-c(1,n)]-coef[2]*temp[-c(n,n-1)]
#temp2new<-temp2[-c(1,2)]-coef[1]*temp2[-c(1,n)]-coef[2]*temp2[-c(n,n-1)]
#trendnew<-trend[-c(1,2)]-coef[1]*trend[-c(1,n)]-coef[2]*trend[-c(n,n-1)]
#partnew<-part[-c(1,2)]-coef[1]*part[-c(1,n)]-coef[2]*part[-c(n,n-1)]
#fitnew = lm(cmortnew~ trendnew + tempnew + temp2new + partnew, na.action=NULL)

#plot(fitnew$residuals)
#acf(fitnew$residuals)
#pacf(fitnew$residuals)



###################################################################
#                          Ejemplo Fuller                         #
###################################################################


###############################
#Se establece el directorio de
#trabajo y se carga la información 
#del indice de producción data[,2]
#y del indice de empleo data[,3]
###############################

setwd("C:/Users/Mayela/Desktop/Expo_series")
data=matrix(scan("2.txt"),ncol=3,byrow=T)


reg<-lm(data[,3]~data[,2])
z<-reg$coefficients[1] +reg$coefficients[2]*data[,2] 
#t<z-seq(19)

plot.ts(data[,3],main="Serie de tiempo",xlab="tiempo",ylab="Indice de empleo",col="blue",xaxt="n")
plot(data[,2],data[,3],xlab="Índice de producción", ylab="Índice de empleo",col="blue",main="Índice de producción vs Índice de empleo")
abline(reg$coefficients,col="red")

axis(1,seq(1,38,4),1955:1964.5)

plot(reg$residuals,main="Residuales",xlab="tiempo",ylab="residuales",col="blue",xaxt="n",type="l")
axis(1,seq(1,38,4),1955:1964.5)

acf(z)
pacf(z)


fit2 <- arima(data[,3], order=c(2,0,0),xreg=cbind(1:length(data[,2]),data[,2]))
datanew<-data