# Script to compute and plot household-wealth-to-gdp ratios
# Espen Henriksen, December 2010

setwd("/Users/espenhenriksen/documents/dropbox/research/global_economy/misc_data_graphs/fact_fiction_figures/")

rm(list=ls())

# Download Table 1.1.5 Gross Domestic Product from BEA and assign it to a data frame called "beatbl115"
beatbl115 <- read.csv("http://www.bea.gov/national/nipaweb/csv/NIPATable.csv?FirstYear=1947&TableName=5&LastYear=2011&ViewSeries=NO&freq=Qtr&3Place=N&Request3Place=N",skip=5,sep=",")
# Name rows and delete empty row(s)
dimnames(beatbl115)[[1]] <- paste("L", beatbl115[,1],": ", sub("^ +", "", beatbl115[,2]), sep = "")
# Delete redudant columns 1 and 2
beatbl115[,2] <- NULL
beatbl115[,1] <- NULL
# Transpose the data set
beatbl115 <- t(beatbl115)
# Format the data set as a time-series data set with quarterly data starting in 1947
beatbl115 <- ts(beatbl115, start=1947, frequency=4)

###

# Download Table 2.1 Personal Income and Its Disposition
beatbl21 <- read.csv("http://www.bea.gov/national/nipaweb/csv/NIPATable.csv?FirstYear=1947&TableName=58&LastYear=2011&ViewSeries=NO&freq=Qtr&3Place=N&Request3Place=N",skip=5,sep=",",na.strings=" --- ")
# Name rows and delete empty row(s)
dimnames(beatbl21)[[1]] <- paste("L",  gsub("(^ +)|( +$)", "", beatbl21[,1]),": ", sub("^ +", "", beatbl21[,2]), sep = "")
beatbl21 <- beatbl21[-c(35,37,39,43,46,47,48,49,50),]
# Delete redudant columns 1 and 2
beatbl21[,2] <- NULL
beatbl21[,1] <- NULL
# Transpose the data set
beatbl21 <- t(beatbl21)
# Format the data set as a time-series data set with quarterly data starting in 1947
beatbl21 <- ts(beatbl21, start=1947, frequency=4)

###

# Table 3.1. Government Current Receipts and Expenditures 
beatbl31 <- read.csv("http://www.bea.gov/national/nipaweb/csv/NIPATable.csv?FirstYear=1947&TableName=86&LastYear=2011&ViewSeries=NO&freq=Qtr&3Place=N&Request3Place=N",skip=5,sep=",",na.strings=" --- ")
# Name rows and delete empty row(s)
dimnames(beatbl31)[[1]] <- paste("L",  gsub("(^ +)|( +$)", "", beatbl31[,1]),": ", sub("^ +", "", beatbl31[,2]), sep = "")
beatbl31 <- beatbl31[-c(30),]
# Delete redudant columns 1 and 2
beatbl31[,2] <- NULL
beatbl31[,1] <- NULL
# Transpose the data set
beatbl31 <- t(beatbl31)
# Format the data set as a time-series data set with quarterly data starting in 1947
beatbl31 <- ts(beatbl31, start=1947, frequency=4)

###

# Download and format Table 5.1. 
beatbl51 <- read.csv("http://www.bea.gov/national/nipaweb/csv/NIPATable.csv?FirstYear=1947&TableName=137&LastYear=2011&ViewSeries=NO&freq=Qtr&3Place=N&Request3Place=N",skip=5,sep=",",na.strings=" --- ")
# Name rows and delete empty row(s)
dimnames(beatbl51)[[1]] <- paste("L",  gsub("(^ +)|( +$)", "", beatbl51[,1]),": ", sub("^ +", "", beatbl51[,2]), sep = "")
beatbl51 <- beatbl51[-c(44,67,68,69,70,71),]
# Delete redudant columns 1 and 2
beatbl51[,2] <- NULL
beatbl51[,1] <- NULL
# Transpose the data set
beatbl51 <- t(beatbl51)
# Format the data set as a time-series data set with quarterly data starting in 1947
beatbl51 <- ts(beatbl51, start=1947, frequency=4)


# Download "Table B.100 Balance Sheet of Households and Nonprofit Organizations" from the Federal Reserve Flow of Funds webpage
download.file(url="http://www.federalreserve.gov/releases/z1/current/Disk/btabs.zip",destfile="btabs.zip",method="auto",quiet=TRUE)
btab100d <- read.csv(unzip(zipfile = "btabs.zip",files="btab100d.prn",list = FALSE),header=TRUE,sep=" ")
btab100d <- ts(btab100d,start=1952, frequency=4)

###

# FL152090005.Q is Line 42 "Net worth"
# FL155035005.Q is Line 3  "Real estate assets"
# FL155012605.Q is Line 43 "Replacement-cost value of structures: Residential"
# FL155012605.Q is Line 46 "Replacement-cost value of structures: Nonresidential (nonprofit)"
# FL153064105.Q is Line 24 "Financial assets: Corporate equities"
# FL153064105.Q is Line 25 "Financial assets: Mutual fund shares"

nwy  <-  btab100d[,"FL152090005.Q"]/beatbl115[,1]/1000
nwc  <-  btab100d[,"FL152090005.Q"]/beatbl115[,2]/1000
hmvy <-  btab100d[,"FL155035005.Q"]/beatbl115[,1]/1000
hrcy <- (btab100d[,"FL155012605.Q"]+btab100d[,"FL165013665.Q"])/beatbl115[,1]/1000
eqy  <- (btab100d[,"FL153064105.Q"]+btab100d[,"FL153064205.Q"])/beatbl115[,1]/1000
tsr   <- (beatbl115[,1]-beatbl115[,2]-beatbl115[,21])/beatbl115[,1]
psr   <- ((beatbl115[,1]-beatbl115[,21])-beatbl115[,2])/(beatbl115[,1]-beatbl115[,21])

# Plot households' net worth over GDP, incl. main components
nwplot <- cbind(nwy,hmvy,hrcy,eqy)
dimnames(nwplot)[[2]] <- c("Households' net worth over GDP","HHs' real-estate assets over GDP","HHs' replacement-cost value of structures over GDP","HHs' corporate equities plus mutual fund shares over GDP")
ts.plot(nwplot,ylim=c(0,5.0),gpars=list(xlab="",ylab="Ratio to GDP",lwd=2.0,lty=c(1,2,2,2),col=c(1:4),axes=F))
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources: Bureau of Economic Analysis, Federal Reserve Flow of Funds Accounts",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(nwplot),cex=0.90,lwd=2.0,lty=c(1,2,2,2),col=c(1:4),bty="n")
dev.print(device=postscript, "nwgdpratio.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "nwgdpratio.pdf")

# Plot the ratio of the worth to GDP and to private consumption, respectively
nwplot <- cbind(nwc,nwy)
dimnames(nwplot)[[2]] <- c("Ratio to private consumption","Ratio to GDP")
ts.plot(nwplot,ylim=c(3.0,7.0),gpars=list(xlab="",ylab="Households' net worth",lwd=2.0,lty=c(1,1),col=c("forestgreen",1),axes=F))
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources: Bureau of Economic Analysis, Federal Reserve Flow of Funds Accounts",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(nwplot),cex=.9,lwd=2.0,lty=c(1,1),col=c("forestgreen",1),bty="n")
dev.print(device=postscript, "nwcnwyratios.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "nwcnwyratios.pdf")

# Plot savings rates
nwplot <- cbind(psr,tsr)
dimnames(nwplot)[[2]] <- c("Private","Total")
ts.plot(nwplot,gpars=list(xlab="",ylab="Savings rate",lwd=2.0,lty=c(1:4),col=c(1:4)))
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources: Bureau of Economic Analysis, Federal Reserve Flow of Funds Accounts",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(nwplot),cex=.9,lwd=2.0,lty=c(1:4),col=c(1:4))
dev.print(device=postscript, "savingsrates.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "savingsrates.pdf")

# Plot consumption and investment ratios
pcty <- beatbl115[,2]/beatbl115[,1]
pity <- beatbl115[,7]/beatbl115[,1]
gcty <- beatbl31[,16]/beatbl115[,1]
gity <- beatbl31[,35]/beatbl115[,1] 
nwplot <- cbind(pcty,pity,gcty,gity)
dimnames(nwplot)[[2]] <- c("Private consumption","Private investment","Gov't consumption","Gov't investment")
ts.plot(nwplot,gpars=list(xlab="",ylab="Ratio to GDP",lwd=2.0,col=c(4,4,2,2),lty=c(1,2,1,2),axes=F))
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources: Bureau of Economic Analysis",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(nwplot),cex=.9,lwd=2.0,col=c(4,4,2,2),lty=c(1,2,1,2),ncol=2,bty="n")
dev.print(device=postscript, "consandinvrates.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "consandinvrates.pdf")
dev.print(device=png,bg = "transparent", width = 600, filename="consandinvrates.png")


# Plot private and government savings
ps <- beatbl51[,44]/beatbl115[,1]
gs <- beatbl51[,47]/beatbl115[,1]
nwplot <- cbind(ps,gs)
dimnames(nwplot)[[2]] <- c("Gross private savings","Gross government savings")
ts.plot(nwplot,ylim=c(-.07,0.23),gpars=list(xlab="",ylab="Ratio to GDP",lwd=2.0,col=c(4,2),lty=c(1,1),axes=F))
abline(h=0,lty=1,col="gray")
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources: Bureau of Economic Analysis",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(nwplot),cex=.9,lwd=2.0,col=c(4,2),lty=c(1,1),ncol=2,bty="n")
dev.print(device=postscript, "savingsrates.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "savingsrates.pdf")

# Real estate equity
reeq <- 1 - btab100d[,"FL153165105.Q"]/btab100d[,"FL155035015.Q"]
#dimnames(nwplot)[[2]] <- c("Gross private savings","Gross government savings")
ts.plot(reeq,gpars=list(xlab="",ylab="Ratio of households real estate equity to real estate value",lwd=2.0,col=c(4,2),lty=c(1,1),axes=F))
#abline(h=0,lty=1,col="gray")
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources:  Federal Reserve Flow of Funds Accounts",side=4,line=.4,cex=.6,adj=0)
#legend("topleft",legend = colnames(nwplot),cex=.9,lwd=2.0,col=c(4,2),lty=c(1,1),ncol=2,bty="n")
dev.print(device=postscript, "realestateequity.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "realestateequity.pdf")

# Plot private and government savings
hrey  <- btab100d[,"FL155035015.Q"]/beatbl115[,1]/1000
hmy <- btab100d[,"FL153165105.Q"]/beatbl115[,1]/1000
remplot <- cbind(hrey,hmy)
dimnames(remplot)[[2]] <- c("Ratio of housholds real estate value to GDP","Ratio of households mortage debt to GDP")
ts.plot(remplot,ylim=c(0.0,1.75),gpars=list(xlab="",ylab="Ratio to GDP",lwd=2.0,col=c(4,2),lty=c(1,1),axes=F))
#abline(h=0,lty=1,col="gray")
axis(1)
axis(2)
mtext("Updated figure from Backus, Henriksen, Lambert and Telmer \"Current Account: Fact and Fiction\" ",side=4,line=-.2,cex=.6,adj=0)
mtext("Data sources:  Federal Reserve Flow of Funds Accounts",side=4,line=.4,cex=.6,adj=0)
legend("topleft",legend = colnames(remplot),cex=.9,lwd=2.0,col=c(4,2),lty=c(1,1),ncol=1,bty="n")
dev.print(device=postscript, "remortage.eps", onefile=FALSE, horizontal=FALSE)
dev.print(device=pdf, "remortage.pdf")