$title  gravity.GMS

* ------------------------------------------------------
* Per Capita Income and the Demand for Skills
* Journal of International Economics

* Justin Caron, Thibault Fally and James Markusen

* December 2019
* ------------------------------------------------------

* ESTIMATE GRAVITY EQUATIONS


**** RUNNING THE FILE MAY REQUIRE PAUSING DROPBOX

* This file is capable of solving for all sectors at once for OLS
* for Poisson ML, solve per sector (in a loop)

$if not set datadir $set datadir "data\"
$setglobal datadir %datadir%

* -----------------------------------------------------------------
* DEFINE GLOBAL VARIABLES

* select dataset to use
* GTAP5
$if not set ds $set ds gtap5
$include loaddata_gtap5.gms


* trade SPECIFICATION ?
* choices : noatc, atc, ctryatc
$if not set tradespec $set tradespec atc

* ---------------------------------------------------------------
* -- SELECT SECTORS AND REGIONS 
* ---------------------------------------------------------------

set i_(i), r_(r);
i_(i) = no; r_(r) = no;

* using all sectors :
i_(i) = yes;

* using all countries which have distance data
r_(r)$(sum(s,importdist(r,s,"ldist"))) = yes;



$ontext
set manuf(i) / PDR, WHT, GRO, V_F, OSD, C_B,
PFB, OCR, CTL, OAP, RMK, WOL, FRS, FSH, COA, OIL, OMN, CMT, OMT, VOL, MIL, PCR, SGR, OFD,
B_T, TEX, WAP, LEA, LUM, PPP, P_C, CRP, NMM, I_S, NFM, FMP, MVH, OTN, ELE, OME, OMF/;

set     serv(i) service sectors / CMN, DWE,ISR,OBS, OFI,OSG,ROS,WTR,TRD,CNS, OTP, ATP, WTP, GDT, ELY/
        tradables(i) the tradable sectors;

tradables(i) = yes;
tradables(serv) = no;

*i_(manuf) = yes;

$offtext

parameter nbr, nbi;
nbr = card(r_); nbi = card(i_);

alias(r_,s_); alias(i_,j_);

display nbr, nbi;
display r_;

parameter	coeffs parameter to store all sector specific stats
 		phiest
		logphiest
		IM
		EX
		cst;


parameter TcEX, TcExctry;


parameter	trade(i,r,s) value of trade and absorbtion parameter in data (from to) INCLUDING INTERMEDIATE USE;

trade(i_,r_,s_) = btrade(i_,r_,s_);
trade(i_,r_,s_)$sameas(r_,s_) = sum(g,vdfm(i_,g,r_));

parameter logtrade;
logtrade(i_,r_,s_)$trade(i_,r_,s_) = log(trade(i_,r_,s_));


parameter linder;

variable sse_, loglikelihood_;


variables
delta_(*,*), IM_, EX_,
                Logfittrade(i,r,s)   defined as trade from r to s,
			tcost0(i,r,s)
                        TcEx_ asymetric trade cost on exporter side
                        TcExctry_ asymetric trade cost on exporter side (ctry-specific)
			cst_(i)
			linder_(i)
			;

Equations	obj, obj_ml, tcostdef, demand,  budget, tradedef, indexdef, fmc, incomedef, cost, cost_nointer, Absorb,
		absorb_nointer, intercostdef, Cost_leontief , intercostdef_leontief, phidef ;

* NLLS : minimizing LOG error terms on trade only 
obj.. sse_ =e=      sum((i_,r_,s_)$logtrade(i_,r_,s_),  sqr( logtrade(i_,r_,s_) - logfittrade(i_,r_,s_)));

* MAXIMUM LIKELIHOOD - assuming poisson distributions

* poisson log-LIKELIHOOD function

obj_ml .. loglikelihood_ =e=  10e-3 * sum((i_,r_,s_),  trade(i_,r_,s_)*logfittrade(i_,r_,s_) - exp(logfittrade(i_,r_,s_)));


tcostdef(i_,r_,s_) .. tcost0(i_,r_,s_)  =e= 
		 -importdist(r_,s_,"ldist")*delta_(i_,"ldist")    
                +importdist(r_,s_,"contig")* delta_(i_,"contig")    
                +importdist(r_,s_,"comlang_off")* delta_(i_,"comlang_off") 
                +importdist(r_,s_,"colony")* delta_(i_,"colony")    
                +importdist(r_,s_,"homebias")* delta_(i_,"homebias")
                +importdist(r_,s_,"bothsea")* delta_(i_,"bothsea")
                +importdist(r_,s_,"rta")* delta_(i_,"rta")
                +importdist(r_,s_,"comcur")* delta_(i_,"comcur")
                +importdist(r_,s_,"common_legor")* delta_(i_,"common_legor")
               + TcEx_(i_,r_)$(not sameas(r_,s_))
                + TcExctry_(r_)$(not sameas(r_,s_));


Tradedef(i_,r_,s_) ..    

	logfittrade(i_,r_,s_)
	=e= IM_(i_,s_) + EX_(i_,r_) + tcost0(i_,r_,s_) + cst_(i_) 
*	+ linder_(i_) * sqr(log(pcexp(s_)) - log(pcexp(r_))) 
;


delta_.L(i_,"ldist") = 1;
delta_.L(i_,"contig") = 0 ;
delta_.L(i_,"comlang_off") = 0 ;
delta_.L(i_,"colony") = 0 ;
delta_.L(i_,"homebias") = 6;
delta_.L(i_,"rta") = 2;
*delta_.L(i_,"bothland") = 0;
delta_.L(i_,"bothsea") = 0;



* if using a constant :

IM_.FX(i_,"usa") = 0;
EX_.FX(i_,"usa") = 0;
Logfittrade.L(i_,r_,s_) = logtrade(i_,r_,s_);
linder_.L(i_) =0;

TcEx_.L(i_,r_) = 0;
TcExctry_.L(r_) = 0;


model calibNLP /obj, tcostdef,   tradedef/;

model calibNLP_ml /obj_ml, tcostdef,  tradedef /;


* Use PATH solver :
*option nlp = pathnlp;

* OLS
*solve calibNLP using nlp minimizing sse_;


* select set of sectors to solve :

set i__(i);


* only sectors with trade
i__(i_)$sum((r_,s_), trade(i_,r_,s_)) = yes;

*i__("tex")  = yes;
*i__("oil")  = yes;



display i__;

* -------------- -------------------------------------------------------
* solving each sector seperately :

file title_ /'title.cmd'/; title_.nd=0; title_.nw=0;

$set stime %time%


i_(i) =no;

loop(i$i__(i),

i_(i) =yes;

*       Update the title bar 
putclose title_ '@title Estimating  ',i.tl, '  (',(round(100*(ord(i)-1)/card(i))),' %% complete)  ',
		      'Start time: %stime%, Current time: %time% -- Ctrl-S to pause'/;
		      
execute 'title.cmd';



* define treatment of border effect:
* ALLOW FOR ASYMETRIC TRADE COSTS :
$if "%tradespec%"=="noatc" TcEx_.FX(i_,r_) = 0; TcExctry_.FX(r_) = 0;
$if "%tradespec%"=="atc" delta_.FX(i_,"homebias") = 0; TcExctry_.FX(r_) = 0;

$if "%tradespec%"=="ctryatc" TcEx_.FX(i_,r_) = 0;




solve calibNLP_ml using nlp maximizing loglikelihood_;

coeffs(dist,"coeff",i_) = delta_.L(i_,dist) ; 

IM(i_,r_) = IM_.L(i_,r_) ;

EX(i_,r_) = EX_.L(i_,r_) ;

cst(i_) = cst_.L(i_);

linder(i_) = linder_.L(i_);


TcEX(i_,r_)  = TcEx_.L(i_,r_);
TcExctry(r_) = TcExctry_.L(r_) ;


i_(i) =no;

);

parameter nbobs;

i_(i) =yes;

phiest(i_,r_) = sum(s_, exp(ex_.L(i_,s_) + tcost0.L(i_,s_,r_) + cst_.L(i_) ));

parameter phiestIV;
phiestIV(i_,r_) = sum(s_$(not sameas(r_,s_)), exp(ex_.L(i_,s_) + tcost0.L(i_,s_,r_) + cst_.L(i_) ));


logphiest(i_,r_) = log(phiest(i_,r_));

nbobs =   sum((i_,r_,s_)$logtrade(i_,r_,s_),  1);

parameter tcostest;

tcostest(i_,s_,r_) = exp(tcost0.L(i_,s_,r_));

* check properties of poisson that sum of fitted imports and sum of true imports is equal
 parameter chkpoisson;

chkpoisson("import",r_, i_) = sum(s_, exp(logfittrade.L(i_,s_,r_)) - trade(i_,s_,r_)) ;
chkpoisson("export",r_, i_) = sum(s_, exp(logfittrade.L(i_,r_,s_)) - trade(i_,r_,s_)) ;

display chkpoisson;

parameter fittrade;

fittrade("fitted",i_,r_,s_) = exp(logfittrade.L(i_,r_,s_));
fittrade("true",i_,r_,s_) = trade(i_,r_,s_);
fittrade("dom flow",i_,r_,r_) = 1;


execute_unload 'estimates\gravityestimates_%ds%_%tradespec%.gdx', coeffs, im, ex,  nbobs, logphiest, phiest , cst, linder, tcostest, TcEX, TcExctry, fittrade, phiestIV;