########################################################################
# Time Series Identification                                           #
# Charles B. Moss                                                      #
# University of Florida                                                #
# cbmoss@ufl.edu                                                       #
########################################################################

dta <- read.csv("TimeSeries-LogDiff.csv")

########################################################################
# dta[,1] - Year                                                       #
# dta[,2] - Log change in U.S. agricultural exports                    #
# dta[,3] - Log change in U.S. agricultural imports                    #
# dta[,4] - Log change in total U.S. agricultural trade                #
# dta[,5] - Pseudo data 1 AR(1)                                        #
# dta[,6] - Pseudo data 2 AR(2)                                        #
# dta[,7] - Pseudo data 3 MA(1)                                        #
# dta[,8] - Pseudo data 4 MA(2)                                        #
# dta[,9] - Pseudo data 5 ARMA(1,1)                                    #
########################################################################
# Choose the variable to analyze                                       #
#                                                                      #

nvar <- 5

z  <- cbind(dta[7:75,nvar],dta[6:74,nvar],dta[5:73,nvar],dta[4:72,nvar],
        dta[3:71,nvar],dta[2:70,nvar],dta[1:69,nvar])

#                                                                      #
# Compute autocorrelations                                             #
#                                                                      #

v  <- var(z)
d1 <- sqrt(diag(v))
c1 <- solve(diag(d1))%*%v%*%solve(diag(d1))

#                                                                      #
# Compute first partial autocorrelation                                #
#                                                                      #

v2 <- v[cbind(1,3,4,5,6,7),cbind(1,3,4,5,6,7)] - v[cbind(1,3,4,5,6,7),2]%*%
         solve(v[2,2])%*%v[2,cbind(1,3,4,5,6,7)]
d2 <- sqrt(diag(v2))
c2 <- solve(diag(d2))%*%v2%*%solve(diag(d2))

#                                                                      #
# Compute second partial autocorrelation                               #
#                                                                      #

v3 <- v[cbind(1,4,5,6,7),cbind(1,4,5,6,7)]-v[cbind(1,4,5,6,7),cbind(2,3)]%*%
         solve(v[cbind(2,3),cbind(2,3)])%*%v[cbind(2,3),cbind(1,4,5,6,7)]
d3 <- sqrt(diag(v3))
c3 <- solve(diag(d3))%*%v3%*%solve(diag(d3))

#                                                                      #
# Compute third partial autocorrelation                                #
#                                                                      #

v4 <- v[cbind(1,5,6,7),cbind(1,5,6,7)]-v[cbind(1,5,6,7),cbind(2,3,4)]%*%
         solve(v[cbind(2,3,4),cbind(2,3,4)])%*%v[cbind(2,3,4),cbind(1,5,6,7)]
d4 <- sqrt(diag(v4))
c4 <- solve(diag(d4))%*%v4%*%solve(diag(d4))

#                                                                      #
# Compute fourth partial autocorrelation                               #
#                                                                      #

v5 <- v[cbind(1,6,7),cbind(1,6,7)]-v[cbind(1,6,7),cbind(2,3,4,5)]%*%
          solve(v[cbind(2,3,4,5),cbind(2,3,4,5)])%*%v[cbind(2,3,4,5),cbind(1,6,7)]
d5 <- sqrt(diag(v5))
c5 <- solve(diag(d5))%*%v5%*%solve(diag(d5))

#                                                                      #
# Compute fifth partial autocorrelation                                #
#                                                                      #

v6 <- v[cbind(1,7),cbind(1,7)] - v[cbind(1,7),cbind(2,3,4,5,6)]%*%
          solve(v[cbind(2,3,4,5,6),cbind(2,3,4,5,6)])%*%v[cbind(2,3,4,5,6),cbind(1,7)]
d6 <- sqrt(diag(v6))
c6 <- solve(diag(d6))%*%v6%*%solve(diag(d6))

rslt <- cbind(c1[cbind(2,3,4,5,6),1],rbind(c2[2,1],c3[2,1],c4[2,1],c5[2,1],c6[2,1]))
print(rslt)

res1 <- lm(dta[7:75,nvar] ~ dta[6:74,nvar])
res2 <- lm(dta[7:75,nvar] ~ dta[6:74,nvar] + dta[5:73,nvar])
res3 <- lm(dta[7:75,nvar] ~ dta[6:74,nvar] + dta[5:73,nvar] + dta[4:72,nvar])
res4 <- lm(dta[7:75,nvar] ~ dta[6:74,nvar] + dta[5:73,nvar] + dta[4:72,nvar] + dta[3:71,nvar])

print(summary.lm(res1))
print(summary.lm(res2))
print(summary.lm(res3))
print(summary.lm(res4))

err <- matrix(res1$residuals)
z1  <- cbind(err[10:69],err[9:68],err[8:67],err[7:66],err[6:65],err[5:64],
             err[4:63],err[3:62],err[2:61],err[1:60])

vp  <- var(z1)
dvp <- sqrt(diag(vp))
cp  <- solve(diag(dvp))%*%vp%*%solve(diag(dvp))
print(cbind(vp[,1],cp[,1]))

lb <- cp[2:10]^2/(60-1:9)
lb <- 60*(60+2)*sum(lb)
print(lb)