log using t4f2.txt, text replace
*
* This is a replication program (stata 17.0) for:
*
* Kiyota, Kozo and Yoshinori Kurokawa (2021) "Factor Intensity Reversals Redux: Feenstra Is Right!"
* Table 4 and Figure 2
*
* Please note that the different version of the software may generate slightly different results from the published version.
*
clear all
set more off
set matsize 5000
use Kiyota_Kurokawa_2021_FIR, clear

*********************************************************************************************
* Table 4: Manufacturing only
*********************************************************************************************

forvalues i = 1973/2009 {
preserve
display "********** year == `i' **********"
keep if year == `i'
keep if ind >= 3 & ind <= 15
keep prefecture ind KL
drop if prefecture == 0
reshape wide KL, i(ind) j(prefecture)
quietly: spearman KL1-KL47

matrix M = r(Rho)
matrix N = r(Nobs)
*matrix list M
drop ind KL1-KL47
matrix v = vec(M)
matrix n = vec(N)

svmat v, names(pre_ind_corr)
svmat n, names(nobs)
*
* convert symmetry matrix to lower triangle matrix
*
merge using C
drop id _merge
rename pre_ind_corr1 pre_ind_corr
rename nobs1 nobs
replace pre_ind_corr = pre_ind_corr + temp
replace nobs = nobs + temp
drop if temp == .
drop temp
*
*
*
gen     group = 12
replace group = 11 if pre_ind_corr >= 0
replace group = 10 if pre_ind_corr >= 0.1
replace group =  9 if pre_ind_corr >= 0.2
replace group =  8 if pre_ind_corr >= 0.3
replace group =  7 if pre_ind_corr >= 0.4
replace group =  6 if pre_ind_corr >= 0.5
replace group =  5 if pre_ind_corr >= 0.6
replace group =  4 if pre_ind_corr >= 0.7
replace group =  3 if pre_ind_corr >= 0.8
replace group =  2 if pre_ind_corr >= 0.9
replace group =  1 if pre_ind_corr == 1

gen year = `i'
save pre_ind_corr`i', replace
restore
}

use pre_ind_corr1973, clear
forvalues i = 1974/2009 {
append using pre_ind_corr`i'
}

table year group, m

tabstat pre_ind_corr, by(year) stat(mean sd n) col(stat) nototal long
sort year
by year: ci means pre_ind_corr

forvalues i = 1973/2009 {
erase pre_ind_corr`i'.dta
}
*
* 95% confidence interval (help sg159 stata)
*
gen z = (1/2)*log((1+pre_ind_corr)/(1-pre_ind_corr))
gen z_l = z - 1.96*(1/sqrt(nobs-3))
gen z_u = z + 1.96*(1/sqrt(nobs-3))
gen lb = (exp(2*z_l)-1)/(exp(2*z_l)+1)
gen ub = (exp(2*z_u)-1)/(exp(2*z_u)+1)
drop z z_l z_u

forvalues i = 2005/2005 {
display "bootstrap confidence interval year == `i'"
bootstrap r(mean), reps(500) seed(1234): summarize pre_ind_corr if year == `i', detail
estat bootstrap, all
}

*********************************************************************************************
* Figure 2
*********************************************************************************************

keep if year == 1975 | year == 1985 | year == 1995 | year == 2005 | year == 2009
label var pre_ind_corr "Spearman's rho"

kdensity pre_ind_corr, nograph generate(x fx)
kdensity pre_ind_corr if year == 1975, nograph generate(fx75) at(x)
kdensity pre_ind_corr if year == 1985, nograph generate(fx85) at(x)
kdensity pre_ind_corr if year == 1995, nograph generate(fx95) at(x)
kdensity pre_ind_corr if year == 2005, nograph generate(fx05) at(x)
kdensity pre_ind_corr if year == 2009, nograph generate(fx09) at(x)

label var fx75 "1975"
label var fx85 "1985"
label var fx95 "1995"
label var fx05 "2005"
label var fx09 "2009"

line fx75 x, lp(shortdash) || line fx85 x, lp(dash) || line fx95 x, lp(longdash) || line fx05 x, lp(solid) sort ytitle(Density)

graph export f2.pdf, replace

*********************************************************************************************
* Table 4: Excludes machinery industries	
*********************************************************************************************

use Kiyota_Kurokawa_2021_FIR, clear

forvalues i = 1973/2009 {
preserve
display "********** year == `i' **********"
keep if year == `i'
keep if ind >= 3 & ind <= 15
drop if ind == 11 | ind == 12 | ind == 13 | ind == 14
*replace ind = 11 if ind == 12 | ind == 13 | ind == 14
*rename ind ind_old
*gen     ind = 1 if ind_old == 6 | ind_old == 7
*replace ind = 2 if ind_old == 8 | ind_old == 9 |  ind_old == 10
*replace ind = 3 if ind_old == 11 | ind_old == 12 | ind_old == 13 | ind_old == 14
*replace ind = 4 if ind_old == 3 | ind_old == 4 | ind_old == 5 | ind_old == 15
keep prefecture ind K LH
*collapse (sum) K LH, by(prefecture ind)
gen KL = K/LH
drop if prefecture == 0
drop K LH
reshape wide KL, i(ind) j(prefecture)
quietly: spearman KL1-KL47

matrix M = r(Rho)
matrix N = r(Nobs)
*matrix list M
drop ind KL1-KL47
matrix v = vec(M)
matrix n = vec(N)

svmat v, names(pre_ind_corr)
svmat n, names(nobs)
*
* convert symmetry matrix to lower triangle matrix
*
merge using C
drop id _merge
rename pre_ind_corr1 pre_ind_corr
rename nobs1 nobs
replace pre_ind_corr = pre_ind_corr + temp
replace nobs = nobs + temp
drop if temp == .
drop temp
*
*
*
gen     group = 12
replace group = 11 if pre_ind_corr >= 0
replace group = 10 if pre_ind_corr >= 0.1
replace group =  9 if pre_ind_corr >= 0.2
replace group =  8 if pre_ind_corr >= 0.3
replace group =  7 if pre_ind_corr >= 0.4
replace group =  6 if pre_ind_corr >= 0.5
replace group =  5 if pre_ind_corr >= 0.6
replace group =  4 if pre_ind_corr >= 0.7
replace group =  3 if pre_ind_corr >= 0.8
replace group =  2 if pre_ind_corr >= 0.9
replace group =  1 if pre_ind_corr == 1

gen year = `i'
save pre_ind_corr`i', replace
restore
}

use pre_ind_corr1973, clear
forvalues i = 1974/2009 {
append using pre_ind_corr`i'
}

table year group, m

tabstat pre_ind_corr, by(year) stat(mean sd n) col(stat) nototal long
sort year
by year: ci means pre_ind_corr

forvalues i = 1973/2009 {
erase pre_ind_corr`i'.dta
}

*
* 95% confidence interval (help sg159 stata)
*
gen z = (1/2)*log((1+pre_ind_corr)/(1-pre_ind_corr))
gen z_l = z - 1.96*(1/sqrt(nobs-3))
gen z_u = z + 1.96*(1/sqrt(nobs-3))
gen lb = (exp(2*z_l)-1)/(exp(2*z_l)+1)
gen ub = (exp(2*z_u)-1)/(exp(2*z_u)+1)
drop z z_l z_u

forvalues i = 2005/2005 {
display "bootstrap confidence interval year == `i'"
bootstrap r(mean), reps(500) seed(1234): summarize pre_ind_corr if year == `i', detail
estat bootstrap, all
}

*********************************************************************************************
* Table 4: Includes agriculture and mining
*********************************************************************************************

use Kiyota_Kurokawa_2021_FIR, clear

forvalues i = 1973/2009 {
preserve
display "********** year == `i' **********"
keep if year == `i'
keep prefecture ind KL
drop if prefecture == 0
reshape wide KL, i(ind) j(prefecture)
quietly: spearman KL1-KL47

matrix M = r(Rho)
matrix N = r(Nobs)
*matrix list M
drop ind KL1-KL47
matrix v = vec(M)
matrix n = vec(N)

svmat v, names(pre_ind_corr)
svmat n, names(nobs)
*
* convert symmetry matrix to lower triangle matrix
*
merge using C
drop id _merge
rename pre_ind_corr1 pre_ind_corr
rename nobs1 nobs
replace pre_ind_corr = pre_ind_corr + temp
replace nobs = nobs + temp
drop if temp == .
drop temp
*
*
*
gen     group = 12
replace group = 11 if pre_ind_corr >= 0
replace group = 10 if pre_ind_corr >= 0.1
replace group =  9 if pre_ind_corr >= 0.2
replace group =  8 if pre_ind_corr >= 0.3
replace group =  7 if pre_ind_corr >= 0.4
replace group =  6 if pre_ind_corr >= 0.5
replace group =  5 if pre_ind_corr >= 0.6
replace group =  4 if pre_ind_corr >= 0.7
replace group =  3 if pre_ind_corr >= 0.8
replace group =  2 if pre_ind_corr >= 0.9
replace group =  1 if pre_ind_corr == 1

gen year = `i'
save pre_ind_corr`i'_a, replace
restore
}

use pre_ind_corr1973_a, clear
forvalues i = 1974/2009 {
append using pre_ind_corr`i'_a
}

table year group, m

tabstat pre_ind_corr, by(year) stat(mean sd n) col(stat) nototal long
sort year
by year: ci means pre_ind_corr

forvalues i = 1973/2009 {
erase pre_ind_corr`i'_a.dta
}

*
* 95% confidence interval (help sg159 stata)
*
gen z = (1/2)*log((1+pre_ind_corr)/(1-pre_ind_corr))
gen z_l = z - 1.96*(1/sqrt(nobs-3))
gen z_u = z + 1.96*(1/sqrt(nobs-3))
gen lb = (exp(2*z_l)-1)/(exp(2*z_l)+1)
gen ub = (exp(2*z_u)-1)/(exp(2*z_u)+1)
drop z z_l z_u

forvalues i = 2005/2005 {
display "bootstrap confidence interval year == `i'"
bootstrap r(mean), reps(500) seed(1234): summarize pre_ind_corr if year == `i', detail
estat bootstrap, all
}

*********************************************************************************************
* Table 4: Human capital
*********************************************************************************************

use Kiyota_Kurokawa_2021_FIR, clear

forvalues i = 1973/2009 {
preserve
display "********** year == `i' **********"
keep if year == `i'
keep if ind >= 3 & ind <= 15
gen KwL = K/wL
keep prefecture ind KwL
drop if prefecture == 0
reshape wide KwL, i(ind) j(prefecture)
quietly: spearman KwL1-KwL47

matrix M = r(Rho)
matrix N = r(Nobs)
*matrix list M
drop ind KwL1-KwL47
matrix v = vec(M)
matrix n = vec(N)

svmat v, names(pre_ind_corr)
svmat n, names(nobs)
*
* convert symmetry matrix to lower triangle matrix
*
merge using C
drop id _merge
rename pre_ind_corr1 pre_ind_corr
rename nobs1 nobs
replace pre_ind_corr = pre_ind_corr + temp
replace nobs = nobs + temp
drop if temp == .
drop temp
*
*
*
gen     group = 12
replace group = 11 if pre_ind_corr >= 0
replace group = 10 if pre_ind_corr >= 0.1
replace group =  9 if pre_ind_corr >= 0.2
replace group =  8 if pre_ind_corr >= 0.3
replace group =  7 if pre_ind_corr >= 0.4
replace group =  6 if pre_ind_corr >= 0.5
replace group =  5 if pre_ind_corr >= 0.6
replace group =  4 if pre_ind_corr >= 0.7
replace group =  3 if pre_ind_corr >= 0.8
replace group =  2 if pre_ind_corr >= 0.9
replace group =  1 if pre_ind_corr == 1

gen year = `i'
save pre_ind_corr`i'_KwL, replace
restore
}

use pre_ind_corr1973_KwL, clear
forvalues i = 1974/2009 {
append using pre_ind_corr`i'_KwL
}

table year group, m

tabstat pre_ind_corr, by(year) stat(mean sd n) col(stat) nototal long
sort year
by year: ci means pre_ind_corr

forvalues i = 1973/2009 {
erase pre_ind_corr`i'_KwL.dta
}

*
* 95% confidence interval (help sg159 stata)
*
gen z = (1/2)*log((1+pre_ind_corr)/(1-pre_ind_corr))
gen z_l = z - 1.96*(1/sqrt(nobs-3))
gen z_u = z + 1.96*(1/sqrt(nobs-3))
gen lb = (exp(2*z_l)-1)/(exp(2*z_l)+1)
gen ub = (exp(2*z_u)-1)/(exp(2*z_u)+1)
drop z z_l z_u

forvalues i = 2005/2005 {
display "bootstrap confidence interval year == `i'"
bootstrap r(mean), reps(500) seed(1234): summarize pre_ind_corr if year == `i', detail
estat bootstrap, all
}

*********************************************************************************************
* Table 4: Excludes high & low land price prefectures (Tokyo, Kanagawa, Osaka, and Akita)
*********************************************************************************************

use Kiyota_Kurokawa_2021_FIR, clear

forvalues i = 1973/2009 {
preserve
display "********** year == `i' **********"
keep if year == `i'
keep if ind >= 3 & ind <= 15
drop if prefecture == 5 | prefecture == 13 | prefecture == 14 | prefecture == 27
keep prefecture ind KL
drop if prefecture == 0
reshape wide KL, i(ind) j(prefecture)
quietly: spearman KL2-KL43

matrix M = r(Rho)
matrix N = r(Nobs)
*matrix list M
drop ind KL2-KL43
matrix v = vec(M)
matrix n = vec(N)

svmat v, names(pre_ind_corr)
svmat n, names(nobs)
*
* convert symmetry matrix to lower triangle matrix
*
merge using C
drop id _merge
rename pre_ind_corr1 pre_ind_corr
rename nobs1 nobs
replace pre_ind_corr = pre_ind_corr + temp
replace nobs = nobs + temp
drop if temp == .
drop temp
*
*
*
gen     group = 12
replace group = 11 if pre_ind_corr >= 0
replace group = 10 if pre_ind_corr >= 0.1
replace group =  9 if pre_ind_corr >= 0.2
replace group =  8 if pre_ind_corr >= 0.3
replace group =  7 if pre_ind_corr >= 0.4
replace group =  6 if pre_ind_corr >= 0.5
replace group =  5 if pre_ind_corr >= 0.6
replace group =  4 if pre_ind_corr >= 0.7
replace group =  3 if pre_ind_corr >= 0.8
replace group =  2 if pre_ind_corr >= 0.9
replace group =  1 if pre_ind_corr == 1

gen year = `i'
save pre_ind_corr`i', replace
restore
}

use pre_ind_corr1973, clear
forvalues i = 1974/2009 {
append using pre_ind_corr`i'
}

table year group, m

tabstat pre_ind_corr, by(year) stat(mean sd n) col(stat) nototal long
sort year
by year: ci means pre_ind_corr

forvalues i = 1973/2009 {
erase pre_ind_corr`i'.dta
}

*
* 95% confidence interval (help sg159 stata)
*
gen z = (1/2)*log((1+pre_ind_corr)/(1-pre_ind_corr))
gen z_l = z - 1.96*(1/sqrt(nobs-3))
gen z_u = z + 1.96*(1/sqrt(nobs-3))
gen lb = (exp(2*z_l)-1)/(exp(2*z_l)+1)
gen ub = (exp(2*z_u)-1)/(exp(2*z_u)+1)
drop z z_l z_u

forvalues i = 2005/2005 {
display "bootstrap confidence interval year == `i'"
bootstrap r(mean), reps(500) seed(1234): summarize pre_ind_corr if year == `i', detail
estat bootstrap, all
}

log close
