### CPIdataAnalysis.R
# This file contains the commands used to analyze the CPI series.


# Source library
#libname <- "H:\\Research\\Vector ARFIMA\\"
libname <- "/home/rebecca/Research/VARFIMA/"

# Source files
source(paste(libname, "FIVAR_VARFI.R", sep=""))

#### Various CPI inflation measures
cpidata <- read.table(paste(libname,"CPIdata.txt",sep=""),header=TRUE)

## Goods and services
# Make inflation measures
goods <- 1200*(cpidata$Commodities[2:733]-cpidata$Commodities[1:732])/cpidata$Commodities[1:732]
services <- 1200*(cpidata$Services[2:733]-cpidata$Services[1:732])/cpidata$Services[1:732]

# Reduce goods and services to get rid of missing observations
goods <- goods[109:732]
services <- services[109:732]

# Combine the two
gsvector <- c(goods,services)
gsvectorDM <- c(goods-mean(goods),services-mean(services))

cp.gs <- CrossPeriodogram(gsvector, 2)

# GPH estimates (based on a variety of number of points)
gph.d <- FIVAR.GPH(gsvector,2)
print(gph.d)
gph.d <- FIVAR.GPH(gsvector,2, numPoints=60)
print(gph.d)
gph.d <- FIVAR.GPH(gsvector,2, numPoints=90)
print(gph.d)


# Initial estimates based on univariate time series
goods.est <- FIVAR.MLE(goods, 1, determinant=TRUE, verbose=FALSE, demean=TRUE)
services.est <- FIVAR.MLE(services, 1, determinant=TRUE, verbose=FALSE, demean=TRUE, init=list(d=c(gph.d[2])))

# Combine these to get initial parameter values
A1start <- diag(c(goods.est$A1[1,1], services.est$A1[1,1]))
SigmaStart <- diag(c(goods.est$Sigma[1,1], services.est$Sigma[1,1]))
dStart <- c(goods.est$d[1], services.est$d[1])

# Full model estimates - MLE and Whittle
us.est.gsDM <- FIVAR.MLE(gsvector, 2,determinant=TRUE, verbose=FALSE, demean=TRUE, init=list(A1=A1start, d=dStart, Sigma=SigmaStart), computeHessian=TRUE)
us.est.gsDM
diag(solve(-us.est.gsDM$Hessian))
diag(solve(-us.est.gsDM$originalHessian))
whittle.est.gsDM <- FIVAR.Whittle(gsvector, 2, verbose=FALSE, demean=TRUE, init=us.est.gsDM, computeHessian=TRUE)
whittle.est.gsDM
diag(solve(-whittle.est.gsDM$Hessian))
FIVAR.loglikelihood(gsvector, A1=whittle.est.gsDM$A1, Sigma=whittle.est.gsDM$Sigma, d=whittle.est.gsDM$d, 
	Mean=whittle.est.gsDM$Mean, determinant=TRUE)
sowell.est.gsDM <- FIVAR.Sowell(gsvector, 2, demean=TRUE, determinant=TRUE, verbose=FALSE, init=us.est.gsDM)
sowell.est.gsDM


# VARFI, too
us.VARFIest.gsDM <- VARFI.MLE(gsvector, 2,determinant=TRUE, demean=TRUE, verbose=FALSE, init=list(A1=A1start, d=dStart, Sigma=SigmaStart), computeHessian=TRUE)
us.VARFIest.gsDM
diag(solve(-us.VARFIest.gsDM$Hessian))
sowell.VARFIest.gsDM <- VARFI.Sowell(gsvector, 2,determinant=TRUE, demean=TRUE, verbose=FALSE, init=us.VARFIest.gsDM)
sowell.VARFIest.gsDM
whittle.VARFIest.gsDM <- VARFI.Whittle(gsvector, 2, demean=TRUE, verbose=FALSE, init=us.VARFIest.gsDM, computeHessian=TRUE)
whittle.VARFIest.gsDM
diag(solve(-whittle.VARFIest.gsDM$Hessian))

# Comparing periodograms
T <- length(gsvector)/2
FIVARspec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    FIVARspec[,,i] <- FIVAR.spectrum(2*pi*i/T, us.est.gsDM$A1, us.est.gsDM$Sigma, us.est.gsDM$d)
}
whittleFIVARspec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    whittleFIVARspec[,,i] <- FIVAR.spectrum(2*pi*i/T, whittle.est.gsDM$A1, whittle.est.gsDM$Sigma, whittle.est.gsDM$d)
}
VARFIspec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    VARFIspec[,,i] <- VARFI.spectrum(2*pi*i/T, us.VARFIest.gsDM$A1, us.VARFIest.gsDM$Sigma, us.VARFIest.gsDM$d)
}
whittleVARFIspec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    whittleVARFIspec[,,i] <- VARFI.spectrum(2*pi*i/T, whittle.VARFIest.gsDM$A1, us.VARFIest.gsDM$Sigma, us.VARFIest.gsDM$d)
}
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(FIVARspec[1,2,1:311]))))
title("Cross-Periodogram and Implied FIVAR Spectral Density")
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(whittleFIVARspec[1,2,1:311]))))
title("Cross-Periodogram and Implied FIVAR Spectral Density (Whittle)")
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(VARFIspec[1,2,1:311]))))
title("Cross-Periodogram and Implied VARFI Spectral Density")
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(whittleVARFIspec[1,2,1:311]))))
title("Cross-Periodogram and Implied VARFI Spectral Density (Whittle)")

# VARFI's linear combination that is white noise
wn <- (goods[2:624]-mean(goods))-.3027*(goods[1:623]-mean(goods))-.4245*(services[1:623]-mean(services))
ts.plot(wn)
title("Implied White Noise Series")
cp.wn <- CrossPeriodogram(wn, 1)
plot(log(Re(cp.wn)))
title("Logged Periodogram of Implied White Noise Series")
acf(wn)



# Checking Sowell computing time
system.time(slik <- VARFI.SowellLogLikelihood(gsvectorDM, us.VARFIest.gsDM$A1, us.VARFIest.gsDM$Sigma, us.VARFIest.gsDM$d))
system.time(ourlik <- VARFI.loglikelihood(gsvectorDM, us.VARFIest.gsDM$A1, us.VARFIest.gsDM$Sigma, us.VARFIest.gsDM$d))
system.time(slikFIVAR <- FIVAR.SowellLogLikelihood(gsvectorDM, us.est.gsDM$A1, us.est.gsDM$Sigma, us.est.gsDM$d))
system.time(ourlikFIVAR <- FIVAR.loglikelihood(gsvectorDM, us.est.gsDM$A1, us.est.gsDM$Sigma, us.est.gsDM$d))



# Fit a VAR instead
library(vars)
gsmatrix <- matrix(nrow=length(goods), ncol=2)
gsmatrix[,1] <- goods-mean(goods)
gsmatrix[,2] <- services-mean(services)
colnames(gsmatrix) <- c("Goods", "Services")
var1 <- VAR(gsmatrix,p=1,type="const")
var2 <- VAR(gsmatrix,p=2,type="const")
# AIC chooses 10 lags
varN <- VAR(gsmatrix,type="const", lag.max=20,ic="AIC")

# Compute the spectral density of the VAR(10)
Alist <- as.list(1:10)
for(i in 1:10){
   nextA <- matrix(nrow=2,ncol=2)
   nextA[1,1:2] <- coef(varN)$Goods[(2*i-1):(2*i)]
   nextA[2,1:2] <- coef(varN)$Services[(2*i-1):(2*i)]
   Alist[[i]] <- nextA
}
sigma2 <- summary(varN)$covres
bigParams <- VARtoVAR1(Alist=Alist, sigma2)

VAR10spec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    VAR10spec[,,i] <- FIVAR.spectrum(2*pi*i/T, bigParams$bigA, bigParams$bigSigma, rep(0,20))[1:2,1:2]
}
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(VAR10spec[1,2,1:311]))))
title("Cross-Periodogram and Implied VAR(10) Spectral Density")

# Compute the covariances/spectral density of the VAR(2)
Alist <- as.list(1:2)
for(i in 1:2){
   nextA <- matrix(nrow=2,ncol=2)
   nextA[1,1:2] <- coef(var2)$Goods[(2*i-1):(2*i)]
   nextA[2,1:2] <- coef(var2)$Services[(2*i-1):(2*i)]
   Alist[[i]] <- nextA
}
sigma2 <- summary(var2)$covres
bigParams <- VARtoVAR1(Alist=Alist, sigma2)

VAR2spec <- array(dim=c(2,2,T))
# (Off by one ordinate from the traditional cross-periodogram)
for(i in 1:T){
    VAR2spec[,,i] <- FIVAR.spectrum(2*pi*i/T, bigParams$bigA, bigParams$bigSigma, rep(0,4))[1:2,1:2]
}
ts.plot(ts(data=log(Mod(cp.gs[2:312,1,2]))), ts(data=log(Mod(VAR2spec[1,2,1:311]))))
title("Cross-Periodogram and Implied VAR(2) Spectral Density")