#  RRR_discussion_Nov_2012.R
#  Pictures for discussion of paper by Pau Rabanal and Juan Rubio-Ramirez,
#  "Low-freq movements of real exchange rates," NYU/AtlFed conference.  
#  See:  http://www.frbatlanta.org/news/conferences/12intl_economics.cfm 
#  Data from FRED, input with functuions written by Paul Backus 
#  Written by:  Dave Backus, November 2012 
# ------------------------------------------------------------------------------
rm(list=ls())

# remember to change \ to /
dir = "C:/Users/dbackus/Documents/Papers/Discussions_and_talks/AtlFed_12"
setwd(dir) 

library("XML")    # for FRED input 
library("xts")    # time series package (also used for FRED)

# 1. FRED functions from Paul 
# ------------------------------------------------------------------------------

# FRED functions from Paul 
callFredAPI <- function(call_string, params) {
  api_key <- "055ba538c874e5974ee22d786f27fdda"  
  url <- paste(
    "http://api.stlouisfed.org/fred/", # base url
    call_string, # subdirectory--documented on fred website
    "?", # separator between web address and parameter list
    paste(
      paste("api_key=", api_key, sep=""),
      paste(
        sapply(
          names(params),
          function(pname) {
            paste(pname, "=", params[[pname]], sep="")
          }
        ),
        collapse="&" # use collapse instead of sep to flatten list
      ),
      sep="&"
    ),
    sep=""
  )
  return(xmlTreeParse(url, useInternal=TRUE))
}

# Helper function to extract specific attributes from FRED's XML
collectAttrs <- function(doc, tag, attr) {
  sapply(
    getNodeSet(doc, paste("//", tag)),
    function(el) { xmlGetAttr(el, attr) }
  )
}

# Downloads the specified series and returns it as a vector, with the dates
# of each observation stored in the vector's names attribute
getFredData <- function(series_id) {
  doc <- callFredAPI("series/observations", list(series_id=series_id))
  dataseries <- as.numeric(collectAttrs(doc, "observation", "value"))
  names(dataseries) <- collectAttrs(doc, "observation", "date")
  return(dataseries)
}

# Downloads the metadata of a FRED series and returns a particular attribute
# A list of available attributes can be found at
#   http://api.stlouisfed.org/docs/fred/series.html
getFredMetadata <- function(series_id, attribute) {
  doc <- callFredAPI("series", list(series_id=series_id))
  attrs <- collectAttrs(doc, "series", attribute)
  return(attrs)
}

# Returns a multivariate time series of the variables specified in series_ids
getFredTable <- function(series_ids) {
  data <- do.call(
    merge,
    lapply(
      series_ids, 
      function(series) { as.xts(getFredData(series)) }
    )
  )
  colnames(data) <- series_ids
  return(data)
}

#  2. Data input 
# ------------------------------------------------------------------------------

# monthly data:  real exchange rate 
fred.sym.m <- c("TWEXMPA")  # real exchange rate 
fred.m <- getFredTable(fred.sym.m)
# kill off some obs
fred.m <- fred.m["1980-01-01/"]  # slash means go to end: "n1/n2", "/n2", "n1/"
fred.m$rer <- log(fred.m$TWEXMPA/100)

#plot(fred.m$rer)

# qrtly data:  net exports 
fred.sym.q <- c("NETEXP", "GDP", "EXPGS", "EXPGSC1", "IMPGS", "IMPGSC1")
fred.q <- getFredTable(fred.sym.q)

# kill off some obs
fred.q <- fred.q["1960-01-01/"]  # slash means go to end: "n1/n2", "/n2", "n1/"
fred.q$nxy <- fred.q$NETEXP/fred.q$GDP
fred.q$tot <- log((fred.q$EXPGSC1/fred.q$EXPGS)/(fred.q$IMPGSC1/fred.q$IMPGS))

#test <- data.frame(fred.q)

#plot(fred.q$nxy, auto.grid="FALSE", minor.ticks="FALSE", major.format="%Y", 
#      main="", ylab="")
#plot(fred.q$tot, auto.grid="FALSE", minor.ticks="FALSE", major.format="%Y", 
#     main="", ylab="")

#  3. Filtering  
# ------------------------------------------------------------------------------

# From Farnsworth, "Econometrics in R," page 25
# http://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf
hpfilter <- function(x,lambda=1600){
  eye <- diag(length(x))
  result <- solve(eye+lambda*crossprod(diff(eye,lag=1,d=2)),x)
  return(result)
}

# monthly data:  real exchange rate  
fred.m$rer.hp <- hpfilter(fred.m$rer) 
fred.m$rer.dev <- fred.m$rer - fred.m$rer.hp 

plot(fred.m$rer,
     auto.grid="FALSE", minor.ticks="FALSE", major.format="%Y", 
     main="", ylab="Log of Real Exchange Rate", 
     mar=c(2,3,2,2), mgp=c(2.5,1,0), col="blue") 
lines(fred.m$rer, col="blue", lwd=2)
lines(fred.m$rer.hp, col="magenta", lwd=2)
abline(h=0, lwd=0.25)

dev.copy2eps(device=postscript, file="RER.eps", width=8, height=6)

# quarterly data:  net exports 
fred.q$nxy.hp  <- hpfilter(fred.q$nxy) 
fred.q$nxy.dev <- fred.q$nxy - fred.q$nxy.hp 

plot(fred.q$nxy,
     auto.grid="FALSE", minor.ticks="FALSE", major.format="%Y", 
     main="", ylab="Ratio of Net Exports to GDP", 
     mar=c(2,3,2,2), mgp=c(2.5,1,0), col="blue") # can't make lwd=2??
lines(fred.q$nxy, col="blue", lwd=2)
lines(fred.q$nxy.hp, col="magenta", lwd=2)
abline(h=0, lwd=0.25)

dev.copy2eps(device=postscript, file="NXY.eps", width=8, height=6)

# terms of trade  
fred.q$tot.hp  <- hpfilter(fred.q$tot) 
fred.q$tot.dev <- fred.q$tot - fred.q$tot.hp 

# cool function from Espen to check column classes (watch for dreaded factors)
frameClasses <- function(x) {unlist(lapply(unclass(x),class))}
frameClasses(fred.q)

fred.ccf <- fred.q #["/1990-01-01"]

# re ccf and xts: http://r.789695.n4.nabble.com/understanding-xts-amp-ccf-td4018370.html
nlag <- 12
xc.raw <- ccf(drop(fred.ccf$nxy), drop(fred.ccf$tot), lag.max=nlag, ylab="Cross-Correlation", xlab="Lag k",main="")
xc.dev <- ccf(drop(fred.ccf$nxy.dev), drop(fred.ccf$tot.dev), lag.max=nlag, ylab="",xlab="Lag k",main="")

plot(x=xc.raw$lag/3600, y=xc.raw$acf, 
     ylab="Cross-Correlation", xlab="Lag k of TOT behind NX", main="", 
     ylim=c(-1,1), type="l", 
     mar=c(2,3,2,2), mgp=c(2.5,1,0), 
     lwd=2, col="blue"
) 
abline(v=0, col="red")
abline(h=0, lwd=0.25)
text(-10, -0.8, pos=4, "Raw data")
lines(x=xc.dev$lag/3600, y=xc.dev$acf, col="magenta", lwd=2) 
text(4, 0.7, pos=4, "HP filtered")

dev.copy2eps(device=postscript, file="XC_nxtot.eps", width=8, height=6)


#  4. Consumption, investment, and GDP 
# ------------------------------------------------------------------------------

fred.sym.q <- c("GDPC96", "GPDIC96", "PCECC96")
fred.q <- getFredTable(fred.sym.q)

# kill off some obs
fred.q <- fred.q["1960-01-01/"]

# take logs 
fred.q$y = log(fred.q$GDPC96) 
fred.q$i = log(fred.q$GPDIC96) 
fred.q$c = log(fred.q$PCECC96)

# why do we need this?
fred.q = data.frame(fred.q)

corrxy <- function(x, y, nlag) {
  gx <- diff(x, lag=nlag)
  gy <- diff(y, lag=nlag)
  corrxy <- cor(gx, gy)
}

lags <- c(1, 2, 3, 4, 8, 12)

corr <- lapply(
        lags, 
        function(n) {
          corrxy(fred.q$y, fred.q$c, n) 
        })

barplot(unlist(corr), names.arg=lags,  
        ylab="Correlation", xlab="Differencing Interval in Quarters", 
        col="blue", ylim=c(0,1), 
        mar=c(2,3,2,2), mgp=c(2.5,1,0))

dev.copy2eps(device=postscript, file="corrs_by_lag.eps", width=8, height=6)


#  5. Spectral analysis  
# ------------------------------------------------------------------------------

#spectrum(as.ts(fred.m$rer))

par(mfcol=c(1,1))

spectrum(fred.m$rer, log="no", main=NULL)
spectrum(fred.m$rer.hp, log="no")
spectrum(fred.m$rer.dev, log="no")

xyzzy <- spectrum(fred.q$nxy)