rm(list = ls())

setwd("/Users/espenhenriksen/documents/research/global_economy")

source("gdfuncs.R")

### MAIN PROGRAM ###

# Parametrize model

# Countries and no of countries
countries <- c("CA","DE","FR","IT","JP","US","GB")
nc <- length(countries)

years <- c("1980","1985","1990","1995","2000","2005")
ny <- length(years)


# Preferences
beta  <- 1.10
sigma <- 4.0

# Technology
alpha = 0.33
delta = .125

# Demographics
I0 <-  4
I  <- 20


load("cohorts.Rdata")
load("sprobs.Rdata")
mortality <- sprobs
rm(sprobs)
cohorts <- cohorts[,,7:12]
dimnames(cohorts)[[3]] <- years
mortality <- mortality[,,7:12]
dimnames(mortality)[[3]] <- years


# Solve for market-clearing prices

NFAGDP <- array(0, dim=c(ny, nc), dimnames=list(years, countries))
imfgdp <- extract.imf.gpd(countries,years)
tfp  <- array(0, dim=c(ny, nc), dimnames=list(years, countries))


# benchmark country, the US, #6
bm <- 6

wedge <- 1

for(ctr in 1:2){
  for(ct in 1:length(countries)){
    tmptfp <- tfp
    if(ct != bm){
      lowtfp <- tmptfp
      lowtfp[,ct] <- -4
      NFAtmp <- compute.nfa.positions(alpha,beta,delta,sigma,I0,I,countries,years,cohorts,mortality,lowtfp,wedge)
      flow <- ((NFAtmp[4,ct,]/NFAtmp[4,bm,]) - (imfgdp[,ct]/imfgdp[,bm]))
      print(flow)

      hightfp <- tmptfp
      hightfp[,ct] <- 0.95
      NFAtmp <- compute.nfa.positions(alpha,beta,delta,sigma,I0,I,countries,years,cohorts,mortality,hightfp,wedge)
      fhigh <- ((NFAtmp[4,ct,]/NFAtmp[4,bm,]) - (imfgdp[,ct]/imfgdp[,bm]))
      print(fhigh)

      midtfp <- tmptfp
      midtfp[,ct] = (hightfp[,ct]+lowtfp[,ct])/2
      NFAtmp <- compute.nfa.positions(alpha,beta,delta,sigma,I0,I,countries,years,cohorts,mortality,midtfp,wedge)
      fmid <- ((NFAtmp[4,ct,]/NFAtmp[4,bm,]) - (imfgdp[,ct]/imfgdp[,bm]))
      print(fmid)
    
      diff <- 1
      iter <- 1

      while (diff > 5E-3) {    
        for(yt in 1:length(years)){
          if(fmid[yt]*flow[yt] < 0){
            hightfp[yt,ct] <- midtfp[yt,ct]
            fhigh[yt] <- fmid[yt]
          } else {
            lowtfp[yt,ct] <- midtfp[yt,ct]
            flow[yt] <- fmid[yt]
          }
        }
        midtfp[,ct] <- (hightfp[,ct]+lowtfp[,ct])/2
        NFAtmp <- compute.nfa.positions(alpha,beta,delta,sigma,I0,I,countries,years,cohorts,mortality,midtfp)
        fmid <- (NFAtmp[4,ct,]/NFAtmp[4,bm,]) - (imfgdp[,ct]/imfgdp[,bm])
        diff <- 0
        for(yt in 1:length(years)){
          diff <- diff + abs(fmid[yt])
        }
        print(diff)

        iter <- iter + 1
        if(iter > 10){
          diff <- 1E-8
          print("iter > 10")
        }      
      }
      tfp[,ct] <- midtfp[,ct]
      print(tfp)
    }
  }
}



NFA <- NFAtmp[3:5,,]
NFAGDP <- array(0, dim=c(ny, nc), dimnames=list(years, countries))
NFAGDP <- t(NFA[3,,])
TAB <- rbind(NFAGDP)
barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(-1.6,1.6))

legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"))
dev.print(device = postscript,"nfaovergdp.eps",onefile = FALSE, horizontal = FALSE)
dev.print(device = pdf,"nfaovergdp_tfptarget.pdf")

barplot(NFA[1,1:7,3]/abs(NFA[1,6,3]),col="navy",ylab="NFA position, normalized by the U.S. position",ylim=c(-1,1))
dev.print(device = postscript,"nfaposition.eps",onefile = FALSE, horizontal = FALSE)
dev.print(device = pdf,"nfaposition_tfptarget.pdf")

tmp <- array(0, dim=c(nc,ny,2),dimnames=list(countries,years,c("Model", "IMF")))
for(i in 1:nc){
  for(j in  1:ny){
    tmp[i,j,1] <- NFA[2,i,j]/NFA[2,6,j]
    tmp[i,j,2] <- imfgdp[j,i]/imfgdp[j,6]
  } 
  #assign(paste(countries[i],sep=""),tmp)
  #print(tmp)
}
relgdpcomp <- tmp
rm(tmp,i,j)

TAB <- rbind(t(relgdpcomp[1,,]),t(relgdpcomp[2,,]),t(relgdpcomp[3,,]),t(relgdpcomp[4,,]),t(relgdpcomp[5,,]),t(relgdpcomp[6,,]),t(relgdpcomp[7,,]))
mp <- barplot(TAB, beside = TRUE, axisnames = FALSE,col=c("deepskyblue","indianred"),legend = c("Model","Data"),ylab="GDP relative to the U.S.",args.legend = list(x = "topleft"))
mtext(1, at = mp, text = c("M", "D"),line = 0, cex = 0.40)
at <- t(sapply(seq(1, nrow(TAB), by = 2), function(x) colMeans(mp[c(x, x+1), ])))
mtext(1, at = at, text = rep(countries, 6),line = 1, cex = 0.50)
mtext(1, at = colMeans(mp), text=years, line = 2)
dev.print(device = postscript,"relgdpcomp_tfptarget.eps",onefile = FALSE, horizontal = FALSE)
dev.print(device = pdf,"relgdpcomp_tfptarget.pdf")


barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="TFP rel. to the United States",ylim=c(0,1.6))
legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
dev.print(device = pdf,"reltfp_tfptarget.pdf")


plot(rownames(NFAGDP),NFAGDP[,1],type="l",lwd=2,ylim = c(-.8,.8),ylab = "NFA (% of GDP)",xlab = "")
for(i in 2:nc){
  lines(rownames(NFAGDP),NFAGDP[,i],lwd=2,col = c(i))
}
legend("topleft",legend = countries,col = c(1:nc),lwd = 2,ncol=2)

tmpNFA <- array(0, dim=c(3, 3, ny), dimnames=list(c("NFA","GDP","NFA/GDP"), c("EU4","USA","JP"),years))
tmpNFA[,1,] <- NFA[,2,]+NFA[,3,]+NFA[,4,]+NFA[,7,]

 TAB <- rbind(exp(tfp))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(-1.6,1.6))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(.5,1.5))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(.5,1.5))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(0,1.5))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(0,1.8))
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="TFP rel. to the United States",ylim=c(0,1.8))
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"))
> ?legend
starting httpd help server ... done
> ?legend
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="TFP rel. to the United States",ylim=c(0,1.8))
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="NFA/GDP",ylim=c(0,1.5))
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="TFP rel. to the United States",ylim=c(0,1.8))
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
> barplot(TAB, beside = TRUE,col=c("grey","blue","orange","darkblue","red","darkgreen"), ylab="TFP rel. to the United States",ylim=c(0,1.6))
> legend("topleft", legend = rownames(TAB),, fill=c("grey","blue","orange","darkblue","red","darkgreen"), bty="n")
> dev.print(device = pdf,"reltfp_tfptarget.pdf")
quartz