rm(list=ls())

#setwd("/Users/espenhenriksen/documents/dropbox/research/global_economy/")
setwd("c:/Documents and Settings/dbackus/My Documents/Userdata/Papers/BCH/data")

library(gtools)
library(gdata)

# structural parameters
delta <- 0.06

countries <- c("AUS","CAN","CHN","CH2","GER","FRA","ITA","JPN","NOR","GBR","USA")

# download data
download.file(url="http://pwt.econ.upenn.edu/php_site/pwt63/pwt63_nov182009version.zip",destfile="pwt63_nov182009version.zip",method="auto",quiet=TRUE)
pwt63 <- read.xls(unzip(zipfile = "pwt63_nov182009version.zip", files="pwt63_w_country_names.xls",list = FALSE), sheet=1, skip=0, header=TRUE)

subset <- {data.frame(country <- pwt63$country,
                      isocode <- pwt63$isocode,
                      year    <- pwt63$year,
                      rgdp    <- pwt63$rgdpl,
                      ki      <- pwt63$ki)}
# rgdp is real investment
# ki is investment share of real gdp

pwt.inv <- NULL
pwt.gdp <- NULL
for(i in 1:length(countries)){
  pwt.inv <- cbind(pwt.inv,subset[subset$isocode == countries[i],]$rgdp*subset[subset$isocode == countries[i],]$ki/100)
  pwt.gdp <- cbind(pwt.gdp,subset[subset$isocode == countries[i],]$rgdp)
}
dimnames(pwt.inv)[[2]] <- countries
pwt.inv <- ts(pwt.inv,start=1950,frequency=1)

first <- mat.or.vec(length(countries),1)
for(i in 1:length(countries)){
  j <- 1
  while(is.na(pwt.inv[j,i])){
    j <- j + 1
  }
  first[i] <- j
}
rm(i,j)

g <- mat.or.vec(length(countries),1)
for(i in 1:length(countries)){
  g[i] <- (pwt.inv[22,i]/pwt.inv[first[i],i])^(1/(22-first[i]))
}

K  <- mat.or.vec(dim(pwt.inv)[1]+1,dim(pwt.inv)[2])*NaN
KY <- mat.or.vec(dim(pwt.inv)[1],dim(pwt.inv)[2])*NaN
for(i in 1:dim(pwt.inv)[2]){
  for(j in first[i]:dim(pwt.inv)[1]){
    if(j == first[i]){
      K[j,i] <- pwt.inv[j,i]/((g[i]-1)+delta)
    }
    K[j+1,i] = (1-delta)*K[j,i] + pwt.inv[j,i]
    if(j <= dim(pwt.inv)[1]){
      KY[j,i] <- K[j,i]/pwt.gdp[j,i]
    }
  }
}
    
KY <- ts(KY,start=1950,frequency=1)

ts.plot(KY,gpars=list(xlab="Year",ylab="Capital-to-output ratio",ylim=c(0,6),col=c(1:ncol(KY)),lwd = 1.75,axes=F))
axis(1, labels=dimnames(KY)[[1]], las=1)
axis(2)
legend("topleft",legend = countries,cex=.8,ncol=5,lwd=2.0,col=c(1:ncol(KY)))
mtext("Data source:Penn World Tables",side=4,line=-.2,cex=.4,adj=0)
dev.print(device=pdf, file="pwtKYratios.pdf")