Using: MATLAB R2014a (8.3.0.532) 64-bit]
…in a working directory containing Ossa’s files.
…run “Set 1” (model/Main_RunMe.m) .
Then, in the same directory,
Using: RStudio 0.98.1049 using R 3.0.1 (64-bit)
…run “Set 2” (model/results/r_main.Rmd) .
Set 1 includes calculations run by Ossa himself (“mycalculations”), which provide the model-setup. From here, I set the tariff shocks with “Shocks = (0:10)/100”, which produces shocks of 0, .01, .02, … , .10 . Next, I create blank MATLAB arrays in which to store the results of the program. Finally, I compute the number of possible clubs, and loop through each club. For each club, I have all of those countries who are ‘In the Club’ each increase their tariffs on those countries who are ‘Out of the Club’, and record the impact this has on each country’s Welfare. Finally, I save the results to .mat files.
Set 2 contains R code written entirely to interpret and graph the outputs of Set 1. The first block of code defines some basic user-friendliness functions and loads some external code from the widely-available CRAN network. Secondly, I load exogenous information: which country-labels were used, how countries are to be weighted, and which tariff-shocks we selected. The third block attempts to reconstruct the original baseline tariffs for each country. The fourth block loads the .mat files from Set 1, and reshapes them into ‘long form’, and deposits them into .csv files for easy use by other software packages. The fifth block selects each country’s optimal tariff (ie, the row which contains the highest ‘Welfare’ value).
The sixth block of code in Set 2 produces .pdf graphs of a countries welfare-response (y axis) to each tariff shock (x-axis). Each page is a different club size: on page 1, clubs are of one country, meaning that there are 7 ways of being in the club and 67 = 42 ways of being out of the club (6 “out countries” per each of the 7 clubs). The seventh block calculates the optimal tariff of each country, for a given club size, by fitting a no-intercept polynomial regression ( (y-1) = 0 + Btariff + (B^2) * tariff ), and then performing simple calculus. The resulting .csv files (“InTheClub.csv” and “OutOfClub.csv”) contain the optimal shock, the welfare which that shock achieves (‘WelfareStar’), as well as, separately, the slope of each regression at a 5% shock level.
The eighth block of code in Set 2 uses the exogenous weights defined in block 2 (above) to calculate the welfare of various aggregate regions: The welfare of a country which is in a club by itself, the weighted average welfare of the 6 countries which are out, and the global weighted average welfare of all 7 countries. The ninth and final code block simply writes a .csv file of Ossa’s tariff data, without aggregating or editing it in any way, so that one might aggregate the industry tariffs (of which there are 33) in a different way (or not at all).
MATLAB File: \model\Main_RunMe.m
%---------------------- % OSSA Model %---------------------- clear all close all clc mycalculations %This program performs a number of frequently needed calculations TRADEs %TRADEs(i,j,s) is the value of trade flowing from country i to country j in industry s; see regions.csv and sectors.csv for the list of regions and industries TARIFFs %TARIFFs(i,j,s) is the tariff applied by country j against industry s imports from country i SIGMA %These are the elasticities of substitution from Table 1 of the paper NX %These are aggregate net exports which are basically zero indicating that this is the purged version of the raw data (as explained in the paper) %------------------------------------------------- %Computing the effects of exogenous tariff changes %------------------------------------------------- clear all close all clc % Shocks = (0:10)/100; Shocks = (25:65)/100; % Shocks = [0.5014162, 0.5330983, 0.5408727, 0.4320262, 0.5654997, 0.5451430, 0.5540263]; L = length(Shocks); C = 7; %seven countries/regions % 1 Brazil % 2 China % 3 EU % 4 India % 5 Japan % 6 ROW % 7 US mycalculations % create a place to dump the results BigResults = zeros(C,L,35,(C-1)); %35 happens to be max { rows( nchoosek(1:7,1:7) ) } BigWhoIsIn = ones(35,C,(C-1)) * -1; % later, we will remove all -1's; % separately, keep all the tarrifs dimTARIFFs = size(TARIFFs); AllTariffs = zeros(dimTARIFFs(1),dimTARIFFs(2),dimTARIFFs(3),L); %for ClubSize = 1 % special edit - for simplicity (when only considering %club sizes of one for ClubSize = 1:(C-1) %clubs can be all or one ClubsOfThisClubSize = nchoosek(1:C,ClubSize); dimClubs = size(ClubsOfThisClubSize); %rows and columns of these clubs qClubs = dimClubs(1); WhoIsIn = zeros(qClubs,C); % declare for use - all out for ClubIndex = 1:qClubs Club = ClubsOfThisClubSize(ClubIndex,:); % Keep track of who is in for n = 1:C % ...for each "i" trading partner... if ismember(n, Club) WhoIsIn(ClubIndex,n) = 1; end end BigWhoIsIn(ClubIndex,:,ClubSize) = WhoIsIn(ClubIndex,:); for Shock = 1:L %for each shock... %have the importer shock the other countries by increasing their %tarrifs by Shocks(n) TARIFFCs=TARIFFs; %reset tarrifs for n = 1:C % ...for each "i" trading partner... if not(ismember(n, Club)) %If they aren't in the climate club, they're getting taxed (rows) for o = 1:C %...for each importer... if ismember(o, Club) %only club-members are doing the taxing (columns) % TARIFFCs(n,o,:)=TARIFFCs(n,o,:)+Shocks(Shock); % % additive method TARIFFCs(n,o,:)= (1+TARIFFCs(n,o,:))*(1+Shocks(Shock))-1; end end end end %welfare calculations NXC=zeros(N,1); %NXC are counterfactual aggregate net exports. I use this to purge the raw data from aggregate trade imbalances as described in the main text LAMBDA=LAMBDABAS; %Select LAMBDABAS if you don't want the lobbying weights, and LAMBDAPOL otherwise [GOVERNMENTWELFAREHAT,WELFAREHAT,WAGEHAT,TRADECs,LOBBYWELFAREHAT,EXPENDITUREHAT]=mycounterfactuals(TARIFFCs,NXC,LAMBDA); BigResults(:,Shock,ClubIndex,ClubSize) = WELFAREHAT; %country-welfare-rows, tarrif-effects (columns), by The Club Itself, ClubSize % save also the tarrifs AllTariffs(:,:,:,Shock) = TARIFFCs; end end end save('results\bigresults.mat', 'BigResults') ; save('results\shocks.mat', 'Shocks') ; save('results\basetariffs.mat', 'TARIFFs') ; save('results\alltariffs.mat', 'AllTariffs') ; save('results\whoisin.mat', 'BigWhoIsIn') ; R Markdown File: \model\results\r_main.Rmd
Ossa Model ( Data for TRICE ) ======================================================== Paul Sztorc `r date()` .Rmd R markdown file. Written in R using version 3.0.1 Made with RStudio 0.98.1049 ```{r PreLoad,echo=FALSE,message=FALSE} rm(list=ls()) Use <- function(package) { if(suppressWarnings(!require(package,character.only=TRUE))) install.packages(package,repos="http://cran.case.edu/") require(package,character.only=TRUE) } Pst <- function(...) paste(...,sep="") setwd("C:/Users/ps583/Documents/GitHub/TRICE/model/results") Use('R.matlab') Use('reshape') Use('ggplot2') ``` This is an R Markdown document. ```{r ExogenousLabelsWeightsShocks} #Label Countries DFlabs <- read.csv("regions.csv") names(DFlabs) <- c("country","importer") Weights <- c(.0279, .1540, .121, .0583, .0540, .3917, .1931) # estimate of gdp weights - exogenous names(Weights) <- DFlabs$importer # Which shocks did we use in Matlab? Shocks <- as.vector(readMat("shocks.mat")$Shocks) Suffix <- Pst("_", round(Shocks[which.min(Shocks)]*100,2), "_", round(Shocks[which.max(Shocks)]*100,2), ".csv") ``` ```{r OriginalTarrifs} # A Completely Optional Step for Looking at the unperturbed Tarif # First, # Import and format tarriff data OriginalTarrifs <- readMat("basetariffs.mat") mOriginalTarrifs <- melt(OriginalTarrifs)[,-5] # lose a useless column names(mOriginalTarrifs) <- c("ExporterGettingTaxed","ImporterApplyingTax","Industry","Tarrif") SidewaysTarrifs <- cast(data=mOriginalTarrifs,formula=ExporterGettingTaxed~ImporterApplyingTax+Industry) AverageTarrifs <- cast(data=mOriginalTarrifs,formula=ExporterGettingTaxed~ImporterApplyingTax,fun.aggregate=mean) # # Optional: # write.csv(mOriginalTarrifs,file="original_tarrifs.csv") # write.csv(SidewaysTarrifs,file="original_tarrifs_stacked_sideways.csv") # write.csv(AverageTarrifs,file="original_tarrifs_avg_by_industry.csv") # average tarrif ResultsOT <- data.frame("country"=1:7, "importer"=DFlabs$importer, "origTarrif"=NA) for(i in 1:7) { TempAT <- as.matrix(AverageTarrifs)[-i,i] TempWeights <- Weights[-i] TempWeights <- TempWeights / sum(TempWeights) TempRes <- TempAT %*% TempWeights ResultsOT$origTarrif[i] <- TempRes } # write.csv(ResultsOT,file="original_tarrifs_avg_by_industry_wgt_by_GDP.csv") ``` ```{r LoadMatlabResults} MatData <- readMat("bigresults.mat") MatClubMembership <- readMat("whoisin.mat") # Shape into useful form mDF <- melt(MatData) mDF <- mDF[mDF$value!=0,-6] # remove stuff which never should have been there names(mDF) <- c("country","shock","club","clubsize","welfare") mDF <- merge(mDF,DFlabs) mDF$cgroup <- paste(mDF$country,mDF$club,sep=".") # which version of the club are we in mDF$clubindex <- paste(mDF$clubsize,mDF$club,sep=":") # Add in the actual shocks - more clear mDF <- merge( mDF, data.frame( "rawshock"=Shocks,"shock"=(1:length(Shocks)) ) ) # Club Membership mCM <- melt(MatClubMembership) mCM <- mCM[mCM$value!=-1,-5] # remove stuff which never should have been there # normal names names(mCM) <- c("club", "country", "clubsize","InTheClub") ClubMembers <- cast(mCM, formula = clubsize + club ~ country, value = "InTheClub") # label the countries names(ClubMembers) <- c("clubsize","club",levels(DFlabs$importer)) ClubMembers$clubindex <- paste(ClubMembers$clubsize,ClubMembers$club,sep=":") # Sanity Check ClubMembers Huge <- merge(ClubMembers,mDF) head(Huge) write.csv(mDF,Pst("r_results/MatlabOutput",Suffix)) write.csv(Huge,Pst("r_results/AnnotatedMatlabOutput",Suffix),row.names=FALSE) ``` ```{r MaxCalculatedNumerically} # Maxes # Maxes <- mDF[mDF$ClubStatus=="In"&mDF$clubsize==2,] # Club only Maxes <- mDF[mDF$clubsize==1,] # Club only MaxesCast <- cast(Maxes,fun.aggregate = max, formula = importer ~ ., value = "welfare" ) # get max row only. # remerge with old data names(MaxesCast) <- c("importer","welfare") MaxesFull <- merge( MaxesCast, mDF[, c("welfare", "cgroup", "rawshock")] ) write.csv(MaxesFull, Pst("r_results/MaxTarrifsFromOssa",Suffix)) ``` ```{r Plots} # Graphical Representation of Trade Data # Fix Suffix (for files later) - must end in '.pdf', of course SuffixPDF <- paste( strsplit(Suffix,".",fixed = TRUE)[[1]][1], ".pdf", sep="") Plots <- vector("list",6) for( i in 1:6 ) { # for each club size # Subset the Data Slice <- mDF[mDF$clubsize==i,] # Build the Plot Plots[[i]] <- ggplot(Slice ,aes(y=welfare,x=rawshock,colour=importer)) + geom_point(size=.5) + geom_line(aes(group=cgroup),alpha=.2) + labs(title=paste("Clubs of size",i)) } # Write to File pdf(file=Pst("AllNations",SuffixPDF)) for( i in 1:6 ) print(Plots[[i]]) dev.off() Plots <- vector("list",6) for( i in 1:6 ) { # for each club size # Subset the Data Slice <- mDF[mDF$clubsize==i,] ClubOnly <- Slice[Slice$welfare >= 1, ] # this happens to be always correct (and graphically what we are interested in, anyway) # Make the Plot Plots[[i]] <- ggplot(ClubOnly ,aes(y=welfare,x=rawshock,colour=importer)) + geom_point(size=.5) + geom_line(aes(group=cgroup),alpha=.2) + stat_smooth(aes(fill=importer), method="lm",formula = y~poly(x,2,raw=TRUE)) + labs(title=paste("Clubs of size",i)) } pdf(file=Pst("ClubNationsOnly",SuffixPDF)) for( i in 1:6 ) print(Plots[[i]]) dev.off() Plots <- vector("list",6) for( i in 1:6 ) { # for each club size # Subset the Data Slice <- mDF[mDF$clubsize==i,] NonClubOnly <- Slice[Slice$welfare <= 1, ] # Make the Plot Plots[[i]] <- ggplot(NonClubOnly ,aes(y=welfare,x=rawshock,colour=importer)) + geom_point(size=.5) + geom_line(aes(group=cgroup),alpha=.2) + stat_smooth(aes(fill=importer), method="lm",formula = y~poly(x,2,raw=TRUE)) + labs(title=paste("Clubs of size",i)) } pdf(file=Pst("r_results/NonClubOnly",SuffixPDF)) for( i in 1:6 ) print(Plots[[i]]) dev.off() ``` ```{r OptimalTarrifs} # Calculate Optimal Tarrif and Slope at 10% # get ready to merge this info ShockDf <- data.frame(shock=1:length(Shocks),rawshock=Shocks) LargeDf <- merge(mDF,ShockDf) for( i in unique( LargeDf$clubsize ) ) { # for each club size # Get the data points Slice <- LargeDf[LargeDf$clubsize==i,] N <- nrow(Slice) # Partition by Club-membership ClubOnly <- Slice[Slice$welfare >= 1, ] NonClubOnly <- Slice[Slice$welfare <= 1, ] # Models - NO INTERCEPT m1 <- lm( I(welfare-1) ~ rawshock:importer+I(rawshock^2):importer + 0, data=ClubOnly) # I (y -1) forces origin to be at 0,0 m2 <- lm( I(welfare-1) ~ rawshock:importer+I(rawshock^2):importer + 0, data=NonClubOnly) ThisRowM1 <- data.frame("ClubSize"=rep(i,7), "importer"=DFlabs$importer, "xBeta"=matrix(coef(m1),ncol=2)[,1], "x2Beta"=matrix(coef(m1),ncol=2)[,2], "df"=summary(m1)$df[2], "r2"= summary(m1)$r.squared ) ThisRowM2 <- data.frame("ClubSize"=rep(i,7), "importer"=DFlabs$importer, "xBeta"=matrix(coef(m2),ncol=2)[,1], "x2Beta"=matrix(coef(m2),ncol=2)[,2], "df"=summary(m2)$df[2], "r2"= summary(m2)$r.squared ) # Create, then append the data: # are we first? FirstRow <- i==unique( LargeDf$clubsize )[1] if(FirstRow) InDF <- ThisRowM1 if(!FirstRow) InDF <- rbind(InDF, ThisRowM1) if(FirstRow) OutDF <- ThisRowM2 if(!FirstRow) OutDF <- rbind(OutDF, ThisRowM2) } # Simple slope calculation - first derivative InDF$SlopeAtTen <- ( InDF$xBeta + (2*InDF$x2Beta*.1)) # where x=.1, what is the slope ? OutDF$SlopeAtTen <- ( OutDF$xBeta + (2*OutDF$x2Beta*.1)) # Basic Calculus-based optimization InDF$OptShock <- -InDF$xBeta / (2*InDF$x2Beta) OutDF$PessShock <- -OutDF$xBeta / (2*OutDF$x2Beta) # Basic Calculus-based optimization InDF$OptShock <- -InDF$xBeta / (2*InDF$x2Beta) OutDF$PessShock <- -OutDF$xBeta / (2*OutDF$x2Beta) # Calculate the Actual Optimized Welfare InDF$WelfareStar <- InDF$xBeta*InDF$OptShock + InDF$x2Beta*InDF$OptShock*InDF$OptShock + 1 # we originally subtracted 1 OutDF$WelfareStar <- OutDF$xBeta*OutDF$PessShock + OutDF$x2Beta*OutDF$PessShock*OutDF$PessShock + 1 # Special Request: Value at 5 % # Calculate the Actual Optimized Welfare InDF$Welfare5pct <- InDF$xBeta*0.05 + InDF$x2Beta*0.05*0.05 + 1 # we originally subtracted 1 OutDF$Welfare5pct <- OutDF$xBeta*0.05 + OutDF$x2Beta*0.05*0.05 + 1 # Dump results write.csv( InDF, file=Pst("r_results/InTheClub",Suffix) ) write.csv( OutDF,file=Pst("r_results/OutTheClub",Suffix) ) ``` ```{r GlobalWelfare} OssaTables <- vector("list",length = length(Shocks) ) for(Shock in Shocks ) { # for each importer # # Global Welfare # # TempDf <- LargeDf[LargeDf$clubsize==1 & LargeDf$rawshock==Shock,] # at 10% shock # at 60% # "In" clubs of one OssaTable <- TempDf[TempDf$country==TempDf$club,c("welfare","importer"),][,c(2,1)] OssaTable <- merge(OssaTable,DFlabs) names(OssaTable) <- c("country","InWelfare",'country') # slightly more complicated subset...produces many results needs aggregation Out <- TempDf[TempDf$country!=TempDf$club,] # Other countries - slightly complex Results <- vector(length=7) for(i in 1:7) { TempWel <- Out[Out$country==i,"welfare"] TempWeights <- Weights[-i] TempWeights <- TempWeights / sum(TempWeights) TempRes <- TempWel %*% TempWeights Results[i] <- TempRes } OssaTable$ElseWelfare <- Results # World - very easy, just weighted average (multiply) Results2 <- vector(length=7) for(i in 1:7) { TempWel <- TempDf[TempDf$country==i,"welfare"] TempRes <- TempWel %*% Weights Results2[i] <- TempRes } OssaTable$WorldWelfare <- Results2 Results3 <- vector(length=7) for(i in 1:7) { TempWel <- Out[Out$country==i,"welfare"] Results3[i] <- median(TempWel) } OssaTable$ElseMedian <- Results3 Results4 <- vector(length=7) for(i in 1:7) { TempWel <- Out[Out$country==i,"welfare"] Results4[i] <- mean(TempWel) } OssaTable$ElseSimpAvg <- Results4 # Add to Databse OssaTable$Weights <- Weights ShockIndex <- (1:length(Shocks))[Shocks==Shock] OssaTables[[ShockIndex]] <- OssaTable } # manipulate data GlobalWelfare <- cast( melt(OssaTables) , L1 + country ~ variable) # merge in the raw shocks names(GlobalWelfare)[1] <- "shock" GlobalWelfare <- merge( GlobalWelfare, data.frame( "rawshock"=Shocks,"shock"=(1:length(Shocks)) ) ) # lose some columns GlobalWelfare <- GlobalWelfare[,c(10,2,9,3,5,6)] write.csv( GlobalWelfare, file=Pst("r_results/OssaTable",Suffix), row.names = FALSE ) ``` ```{r TariffAnalysis} AllTariffs <- melt( readMat('alltariffs.mat') )[,-6] names(AllTariffs) <- c("exporter_taxed","importer_taxing","industry","shock","value") ShockDf <- data.frame(shock=1:length(Shocks),rawshock=Shocks) mAllTariffs <- merge(AllTariffs,ShockDf) write.csv(mAllTariffs,Pst("r_results/AllTariffs",Suffix)) ```