/*
Douglas Hibbs 2008-02-09
This program documents the stata code for Hibbs and Olsson, 
PNAS vol.101 No.10, March 9, 2004 pages 3715-3720,
and it generates a data appendix in xls format 
which along with the stata program and two stata data files
are posted at Hibbs' web site: www.douglas-hibbs.com ("scientific articles")
*/
 
 * PROGRAM FOR TRANSITION DATES ANALYSIS, EQ. 3 AND TABLE 2
clear 
version 9
set mem 200m

cd c:\myfiles-2008\pnas2004_2008-02-09-documentation
use transition6_2008-02-09-documentation.dta

describe, fullnames
sort t_nA 
list plants animals A_tilde  t_nA

* TABLE 2, ESTIMATES FOR TRANSITIONS DATES, EQ. 3 (t_nA)
nl ( t_nA = (1/{gamma})* (log((log(1/A_tilde)) / (log({A_A}/A_tilde))) ) ),  ///
variables(t_nA  A_tilde) initial(gamma 0.0005 A_A 10) iterate(500) nolog

predict t_nA_hat
predict residual_t_nA, residuals
list t_nA t_nA_hat residual_t_nA A_tilde

* FIGURE 1
  graph twoway scatter t_nA t_nA_hat A_tilde, sort(t_nA_hat) ///
legend(label(1 "Observed") label(2 "Fitted") position(1) ///
ring(1) region(style(none))) msymbol(circle i) ///
msize(medium) mcolor(black) ///
ytitle("Transition date") xtitle("Biogeographic endowment") ///
connect(i l) xlabel(0 10 to 100) ylabel(1000 2000 to 8000, ///
angle(0) nogrid)xtick(0 10 to 100) ytick(1000 2000 to 8000) jitter(1) ///
graphregion(color(white))


* PROGRAM FOR LONG-RUN DEVELOPMENT ANALYSES, EQ.5, TABLE 3, FIGS.3-4

clear 
version 9
set mem 200m
cd c:\myfiles-2008\pnas2004_2008-02-09-documentation
* log using hibbs-and-olsson_pnas-2004
use longdevN112_2008-02-09-documentation.dta
describe, fullnames

/* NOTE: institutions variable is supplemented from 
data in Kaufmann et al., re-normed to a 0,100 scale.
See (http://econ.worldbank.org/files/11783_wps2772.pdf) 
*/
 
 * THE TIMELINE SPANS 12 MILLENNIA, t=1 at 10,000 BC to t=12,000 at AD 2000

  
/*  TRANSTION6 ESTIMATES FROM EQ 3., TABLE 2
gen A_A6 = 9.2287678
gen gamma6 = 0.0003059741
 */ 

/* GEOGRAPHY
egen maxclimate = max(climate) 
gen relclimate = (climate)/maxclimate


egen maxaxis = max(axis) 
gen relaxis = axis/maxaxis

gen abslatitude = (abs(latitude))/90
egen maxlatitude = max(abslatitude) 
gen rellatitude = abslatitude/maxlatitude

gen A_tilde_neoeurope = A_tilde
replace A_tilde_neoeurope =  ///
(1.0/3)*(rellatitude+relclimate+relaxis)*100 if neoeurope==1
label var A_tilde_neoeurope "actual neoeurope biogeograpy without imputation of Europe max value=100"

gen geography = (1.0/3)*(rellatitude+relclimate+relaxis)*100

gen lny = log(gdppc97)
gen lngeography = log(geography)
*/

/* NOTE: NeoEuropes (Australia, Canada, NZ and US)were imputed the 
max/highest A_tilde (biogeography) score (100) to form the 
Transition Score variable (T - t_An) as  described in the main text 
of Hibbs&Olsson PNAS 2004.  However the neoEuropes were not assigned 
 max A_tilde (biogegraphy) for purposes of the 
A_tilde on Geography regression in Model 1, Table 3. */

/* TABLE 3, MODEL 1
The result for model 1 generated below is stronger than that 
reported in the PNAS article. There must have been an error in transciption.
*/
reg A_tilde_neoeurope geography if smpl_reg==1

* TABLE 3, MODEL 2
reg Tminust_nA lngeography if smpl_reg==1

* TABLE 3, MODELS 3-6
reg lny Tminust_nA  if smpl_reg==1
reg lny lngeography  if smpl_reg==1
reg lny lngeography Tminust_nA Tminust_nA_sq  if smpl_reg==1
reg lny lngeography Tminust_nA Tminust_nA_sq institutions if smpl_reg==1

* TABLE 3, MODELS 7-8
reg institutions geography if smpl_reg==1
reg institutions Tminust_nA if smpl_reg==1

* TABLE 3, MODELS 9
reg lny institutions if smpl_reg==1