### File: PhillipsExample.R
## Copyright (c) 2008, Rebecca Sela
## All rights reserved.
##
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are met:
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.
##     * Redistributions in binary form must reproduce the above copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
##     * Neither the name of the <organization> nor the
##       names of its contributors may be used to endorse or promote products
##       derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY Rebecca Sela ``AS IS'' AND ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL Rebecca Sela BE LIABLE FOR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##
## This code can be cited by citing the paper:
## "Computationally Efficient Gaussian Maximum Likelihood Methods for Vector
## ARFIMA Models" by Sela and Hurvich.
##
## Please report bugs to Rebecca Sela <rsela@stern.nyu.edu>.
##

# This file contains the code used to compute the FIVAR and VARFI models (and
# some associated quantities) for the annual unemployment rate and
# inflation rate data.


########## Annual Phillips curve data (from Wooldridge)
phillipsann <- read.table("phillips.raw",header=TRUE)
Unemployment <- phillipsann$Unemployment
Inflation <- phillipsann$Inflation

# Graphing the cross-correlation function
ccf(Unemployment, Inflation)

# Combine the data into a single column vector
phillips <- as.matrix(phillipsann[,2:3])
dim(phillips) <- NULL

# GPH estimates
gph.d <- FIVAR.GPH(phillips, 2)
print(gph.d)

# Demeaning
phillipsDM <- as.matrix(phillipsann[,2:3])-colMeans(as.matrix(phillipsann[,2:3]))

### Full FIVAR model estimates - MLE and Whittle
us.est <- FIVAR.MLE(phillipsDM, 2,determinant=TRUE, verbose=FALSE)
sowell.est <- FIVAR.Sowell(phillipsDM, 2,determinant=TRUE, verbose=TRUE, init=us.est)
whittle.est <- FIVAR.Whittle(phillipsDM, 2, verbose=FALSE)
us.est
whittle.est

# Log likelihoods (regression approximations) at the MLE and Whittle estimates
FIVAR.loglikelihood(phillipsDM, A1=us.est$A1, Sigma=us.est$Sigma, d=us.est$d,determinant=T)
FIVAR.loglikelihood(phillipsDM, A1=whittle.est$A1, Sigma=whittle.est$Sigma, d=whittle.est$d,determinant=T)

# Autocovariances implied by the MLE estimates
covsFIVAR <- FIVAR.autoCov(d=us.est$d,Sigma=us.est$Sigma,us.est$A1,n=49)

### Full VARFI model estimates - MLE and Whittle
# This one doesn't converge to the right values
us.VARFIest <- VARFI.MLE(phillipsDM, 2,determinant=TRUE, verbose=FALSE)
# This one does
us.VARFIest2 <- VARFI.MLE(phillipsDM, 2,determinant=TRUE, init=us.est, verbose=FALSE)
sowell.VARFIest2 <- VARFI.Sowell(phillipsDM, 2,determinant=TRUE, init=us.VARFIest2, verbose=FALSE)
whittle.VARFIest <- VARFI.Whittle(phillipsDM, 2, verbose=FALSE)
us.VARFIest
whittle.VARFIest

# Log likelihoods - with the regression approximation and with Sowell's exact algorithm
VARFI.loglikelihood(phillipsDM, A1=us.VARFIest$A1, Sigma=us.VARFIest$Sigma, d=us.VARFIest$d,determinant=T)
VARFI.loglikelihood(phillipsDM, A1=whittle.VARFIest$A1, Sigma=whittle.VARFIest$Sigma, d=whittle.VARFIest$d,determinant=T)
VARFI.SowellLogLikelihood(phillipsDM, A1=us.VARFIest$A1, Sigma=us.VARFIest$Sigma, d=us.VARFIest$d,determinant=T)
VARFI.SowellLogLikelihood(phillipsDM, A1=whittle.VARFIest$A1, Sigma=whittle.VARFIest$Sigma, d=whittle.VARFIest$d,determinant=T)
VARFI.SowellLogLikelihood(phillipsDM, A1=sowell.VARFIest2$A1, d=sowell.VARFIest2$d, Sigma=sowell.VARFIest2$Sigma)

# Autocovariances implied by the MLE estimates
covsVARFI <- VARFI.autoCov(d=us.VARFIest$d,Sigma=us.VARFIest$Sigma,us.VARFIest$A1,n=49)
covest <- VARFI.autoCov(F=us.VARFIest2$A1, d=us.VARFIest2$d, Sigma=us.VARFIest2$Sigma,13)

### Fitting a VAR(2) for comparison
library(vars)
var1 <- VAR(phillipsann[,2:3],p=1,type="const")
var2 <- VAR(phillipsann[,2:3],p=2,type="const")
# AIC chooses 2 lags
varN <- VAR(phillipsann[,2:3],type="const", lag.max=10,ic="AIC")