$title  CF_skill.GMS

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

* Justin Caron, Thibault Fally and James Markusen

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

* COUNTERFACTUAL SIMULATIONS WITH COMIN ET AL. 2016 (CLM) PREFERENCES

*-----------
* homothetic or non-homothetic preferences? choices: NH, H
$if not set spec $set spec NH

* demand specification? choices: tc, theta4
$if not set demandest $set demandest theta4

* --- CHOICE OF SCENARIO:

* to run the UNIFIED counterfactual, set tradetarg and prodtarg to yes (targets gdp and export growth)
* to run the TRADE-ONLY counterfactual, set tradetarg to load and prodtarg to no
* to run the GROWTH-ONLY counterfactual, set prodtarg to load and tradetarg to no

$if not set tradetarg $set tradetarg yes
*$if not set tradetarg $set tradetarg no
*$if not set tradetarg $set tradetarg load


$if not set prodtarg $set prodtarg yes
*$if not set prodtarg $set prodtarg no
*$if not set prodtarg $set prodtarg load



* --- NB OF ITERATIONS FOR GE CONVERGENCE:
Sets
$if "%prodtarg%" == "yes"  T GE iteration  /T1*T6/
$if "%prodtarg%" == "no"  T GE iteration  /T1*T11/
$if "%prodtarg%" == "load"  T GE iteration  /T1*T11/

* --- NB OF ITERATIONS FOR GE CONVERGENCE:

$if "%prodtarg%" == "yes" TARG targetting iterations /targ1*targ16/;
$if "%prodtarg%" == "load"  TARG targetting iterations /targ1*targ1/;
$if "%tradetarg%" == "load" TARG targetting iterations /targ1*targ1/;

* note: we are using 16 convergence for final runs but results are very close with 5-6 iterations

* --- FOR DECOMPOSITION COUNTERFACTUALS, LOAD SHOCKS FROM UNIFIED CF
parameter reporting stores all model results
          reporting_ loaded model results from unified cf;

$if exist 'results/%demandest%_CLM/CF_Skill_prod_yes_trade_yes_%spec%_CLM.gdx'  $gdxin results/%demandest%_CLM/CF_Skill_prod_yes_trade_yes_%spec%_CLM.gdx
$if exist 'results/%demandest%_CLM/CF_Skill_prod_yes_trade_yes_%spec%_CLM.gdx'  $load reporting_ = reporting

* ------------------------------------------------------------------------------------------------------
* ------------------------------------------------------------------------------------------------------

Sets
r country  /
col,aus,
nzl,chn,hkg,jpn,kor,twn,idn,mys,phl,sgp,tha,vnm,bgd,ind,lka,can,usa,mex,
per,ven,arg,bra,chl,ury,aut,bel,dnk,fin,fra,deu,gbr,grc,irl,ita,lux,nld,prt,esp,swe,che,
alb,bgr,hrv,cze,hun,mlt,pol,rom,svk,svn,est,lva,ltu,rus,cyp,tur,mar,bwa,mwi,moz,tza,zmb,zwe,uga/;
alias (r,s);


* DEFINE TARGETS for unified model run
parameters income_target, tcost_target, target_rescale_factor factor by which target is rescaled (to help model solve);
tcost_target(r) = 1;
income_target(r) = 1;
target_rescale_factor = 10;

*  Load the GDP and export/GDP growth rates used for targetting
$include Targets.gms

* trade targets (export/gdp ratios)
* USING PWT9 (Penn tables) for trade/gdp growth
$if "%tradetarg%" == "yes" tcost_target(r)$exp_gdp(r,"1995") = ((exp_gdp(r,"2010") / exp_gdp(r,"1995")-1 )/target_rescale_factor) +1;
$if "%tradetarg%" == "yes" tcost_target("SGP")$exp_gdp("SGP","1995") = ((exp_gdp("SGP","2010") / exp_gdp("SGP","1995")-1 )/(1.5*target_rescale_factor)) +1;

* if using WB
*$if "%tradetarg%" == "yes" tcost_target(r)$exp_gdp_WB(r,"1995") = ((exp_gdp_WB(r,"2010") / exp_gdp_WB(r,"1995")-1 )/target_rescale_factor) +1;
* taking TWN from other dataset
*$if "%tradetarg%" == "yes" tcost_target("TWN")= 1.02;

* productivity targets (per capita GDP growth)
$if "%prodtarg%" == "yes" income_target(r)$pcgdp(r,"1995") = ((pcgdp(r,"2010") / pcgdp(r,"1995")-1 )/target_rescale_factor) +1;
$if "%prodtarg%" == "yes" income_target("mwi")= 1.02;
$if "%prodtarg%" == "yes" income_target("ZMB")$pcgdp("ZMB","1995") = ((pcgdp("ZMB","2010") / pcgdp("ZMB","1995")-1 )/(1.8*target_rescale_factor)) +1;
display income_target;

* -------------

Sets
i industry  /

* gtap8
*isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,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,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi

*gtap 5
pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,"for", fsh,col,oil,gas,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,ely,gdt,wtr,cns,trd,otp,wtp,atp,cmn,ofi,isr,obs,ros,osg/;

alias (i,j);

Sets
f factor /Land,UnskLab,SkLab,Capital,NatlRes/;


* -------------
* Loading data for model calibration:
parameter sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, fittedexp, coeffs;
$gdxin estimates/DATA_gtap5_%demandest%_avg.gdx
$load sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, fittedexp, coeffs




************************
* CLM:

theta(i)=4;

* Loading CLM estimated parameter:
parameter sigma_p, epsilon_diff;
$gdxin 'estimates\estimates_gtap5_CLM.gdx'
$load sigma_p, epsilon_diff


*sigma_p(i) = 0.76;

display sigma_p, epsilon_diff;

parameter epsi epsilon parameter in CLM;
epsi(i) = sigma_p + epsilon_diff(i);

display epsi;


$if "%spec%"=="H"  epsi(i) = sigma_p ;

* avoid extreme values:
theta(i)$(theta(i) > 60)  = 60;

*ensure theta's such that the average equals = 4 and min = 1:
theta(i)$(theta(i) = 0)  = 0;
theta(i)$(theta(i) < 0)  = 0;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));
theta(i)$(theta(i) = 1)  = 1;
theta(i)$(theta(i) < 1)  = 1;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));


$if "%demandest%"=="theta4" theta(i) = 4;

* some statistics:
parameter avgtheta, avgsigma, average_mu;

*avgtheta = sum(i$(coeffs("mu","coeff",i)  and not sigmazero(i) and not skipsector(i)),(theta(i) )) / sum(i$(coeffs("mu","coeff",i)  and not sigmazero(i) and not skipsector(i)),1);
avgtheta = sum(i$(coeffs("mu","coeff",i)),(theta(i) )) / sum(i$(coeffs("mu","coeff",i)  ),1);
avgsigma = sum(i$(coeffs("mu","coeff",i)),(sigma(i) )) / sum(i$(coeffs("mu","coeff",i)  ),1);

parameter check_price_elast;
check_price_elast("estimated mu",i) = coeffs("mu","coeff",i);
check_price_elast("implied mu",i)$theta(i) = (sigma(i)-1)/theta(i);

average_mu = sum(i$(coeffs("mu","coeff",i)),(check_price_elast("implied mu",i) )) / sum(i$(coeffs("mu","coeff",i)  ),1);

display sigma, theta, check_price_elast;
display avgtheta, avgsigma, average_mu;


** REDEFINING ABSORPTION AND PRODUCTION FOR CONSISTENCY (GTAP isn't balanced):
parameters production, absorption;

production(i,r) = sum(s, trade("data",i,r,s));
* INSTEAD OF:
*production(i,r) = vom(i,r);

absorption(j,r) = fd(j,r) + sum(i,intersh(j,i,r) * production(i,r));
* INSTEAD OF:
*absorption(i,r) = sum(j, vdfm(i,j,r)+ vifm(i,j,r));
* THIS GUARANTEES WELL DEFINED GOSH COEFFICIENTS -- but may lead to consistency issues with import flows.



** COMPUTING SHARES FOR INITIAL MODEL CALIBRATION:
* demand shares
* trade shares
* factor shares

parameter
dem_sh(i,r)  share of sector i in total final demand
trade_sh(i,r,s)  share of country r in country s's absorption
prod_sh(i,r,s)  share of country s in country r's sales
fact_sh(f,i,r)   share of sector i in demand for factor f
endow_sh(f,r) share of factor f in total endowment (income)
gosh_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_f(j,r) gosh matrix coefficient (how much j is used by final consumers)
gosh_alt_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_alt_f(j,r) gosh matrix coefficient (how much j is used by final consumers);


* j being used in industry i:
gosh_m(j,i,r)$absorption(j,r) = intersh(j,i,r) * production(i,r) / absorption(j,r);

gosh_f(j,r) = 1 - sum(i,gosh_m(j,i,r));
* normalized to unity when absorption is zero

* Alternative def to approx Gosh coef in closed econ
gosh_alt_f(j,r) = (absorption(j,r) - sum(i,intersh(j,i,r) * production(i,r))) / absorption(j,r) ;
gosh_alt_m(j,i,r)$absorption(j,r) = intersh(j,i,r) * absorption(i,r) / (absorption(j,r)*gosh_alt_f(j,r) + sum(i.local, intersh(j,i,r) * absorption(i,r)) );
gosh_alt_f(j,r) = 1 - sum(i,gosh_alt_m(j,i,r));

display gosh_f;
* Some negative values because of rounding errors, not a problem

dem_sh(i,r) =  fd(i,r) / sum(i.local, fd(i,r));

* defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));

parameter comp_ab1, comp_ab2;
comp_ab1(i,s)$absorption(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
comp_ab2(i,s)$(round(sum(r,trade_sh(i,r,s)) - 1,14)) = absorption(i,s);
display comp_ab1, comp_ab2;

trade_sh(i,s,s)$(sum(r.local,trade("data",i,r,s))=0)  = 1;

* defined as: from r to s:
prod_sh(i,r,s)$sum(s.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(s.local , trade("data",i,r,s));

* adjustment when sum is zero:
prod_sh(i,r,r)$(sum(s.local , trade("data",i,r,s))=0)  = 1;

fact_sh(f,i,r) =  beta(f,i, r) * vom(i,r) / sum(i.local,  beta(f,i,r) * vom(i,r))  ;

endow_sh(f,r) = sum(i,vfm(f,i,r)) / sum((f.local,i),vfm(f,i,r)) ;

* ------------------------------------------------------------------------
* CHECKS TO ENSURE THAT CHG=1 IS A SOLUTION:

parameters CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
CHECK_budget(r)= round(sum(i,dem_sh(i,r)) - 1,14);
CHECK_absorption(j,r) = round(gosh_f(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r)) - 1,14);
CHECK_production(i,r) = round(sum(s,prod_sh(i,r,s)) - 1,14);
CHECK_supplier(i,r) = round(prod(f$beta(f,i,r),1**beta(f,i,r))**theta(i) * prod(j,1**intersh(j,i,r))**theta(i) - 1,14);
CHECK_index(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
CHECK_fmc(f,r) = round(sum(i,fact_sh(f,i,r)) - 1,14);
CHECK_income(r) = round(sum(f,endow_sh(f,r)) - 1,14);

display CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;


* -------------------------------------------------------------------------
* CODE TO COMPUTE FITTED PRODUCTION SHARES -- CAN BE SKIPPED ONCE 'fitted_prod_shares_%demandest%.gdx HAS BEEN CREATED
*$ontext


* recreating absorption and production given fitted (final) demand and observed import shares:

* quick model to solve for production, based on fitted consumption and observed trade shares and IO links:
positive variables reconstruct_absorption(i,r) fitted absorption
                   reconstruct_production(i,r) fitted production;

parameter final_demand, trade_sh_adj;
final_demand(i,r) = fittedexp("true",i,r);
trade_sh_adj(i,r,s) = trade("data",i,r,s)/absorption(i,s);


equations eq_IOreconstruct, eq_prodreconstruct;

eq_prodreconstruct(i,r)..   reconstruct_production(i,r) =e= sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption(i,s) );
eq_IOreconstruct(i,r)..     reconstruct_absorption(i,r) =e= final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production(j,r));


model Reconstruct / eq_prodreconstruct.reconstruct_production, eq_IOreconstruct.reconstruct_absorption /;
Reconstruct.iterlim = 1000000;

reconstruct_production.L(i,r) = sum(s.local , trade("data",i,r,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

parameters CHECK_production_bis;

CHECK_production_bis(i,r)$reconstruct_production.L(i,r) = round(1 - sum(s.local , trade("data",i,r,s)) / reconstruct_production.L(i,r), 14);
display CHECK_production_bis;
* not exactly zero but very close

parameters fittedprod, fitted_prod_sh;

fittedprod(i,r,"true") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"true")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / reconstruct_production.L(i,r);

* now: computing how nh pref distord export and production patterns:

* note: problems when fitted final demand is nonzero but true final demand is zero:
final_demand(i,r)$fittedexp("true",i,r) = fittedexp("non-homoth",i,r);

*display reconstruct_production.L, reconstruct_absorption.L;

reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

fittedprod(i,r,"non-homoth") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"non-homoth")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / reconstruct_production.L(i,r);

parameter distortion_prod;
distortion_prod(i,r,"non-homoth")$fittedprod(i,r,"true") = fittedprod(i,r,"non-homoth") / fittedprod(i,r,"true");


* note: problems when fitted final demand is nonzero but true final demand is zero:
final_demand(i,r)$fittedexp("true",i,r) = fittedexp("homoth",i,r);

reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

fittedprod(i,r,"homoth") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"homoth")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / reconstruct_production.L(i,r);

distortion_prod(i,r,"homoth")$fittedprod(i,r,"true") = fittedprod(i,r,"homoth") / fittedprod(i,r,"true");


display distortion_prod;

execute_unload 'fitted_prod_shares_%demandest%.gdx', fittedprod, fitted_prod_sh, distortion_prod;

*$exit

*$offtext


*$gdxin fitted_prod_shares_%demandest%.gdx
*parameters fittedprod, fitted_prod_sh;

*$load fitted_prod_sh



***
*** SPECIFICATION:

$if "%spec%"=="H" prod_sh(i,r,s) = fitted_prod_sh(i,r,s,"homoth");
$if "%spec%"=="NH" prod_sh(i,r,s) = fitted_prod_sh(i,r,s,"non-homoth");
***

$if "%spec%"=="NH" dem_sh(i,r) = fittedexp("non-homoth",i,r) / sum(i.local, fittedexp("non-homoth",i,r));
$if "%spec%"=="H" dem_sh(i,r) = fittedexp("homoth",i,r) / sum(i.local, fittedexp("homoth",i,r));

* note: employment shares and income shares aren't adjusted as they are less likely to matter


***
* ALLOWING FOR CES PRODUCTION FUNCTION BETWEEN LABOR TYPES:

Set f_noLab factors other than labor /Land,Capital,NatlRes/;

Set f_Lab types of labor /SkLab,UnskLab/;


parameter lab_sh(f_Lab,i,r), betalab, rho elasticity of substitution between labor types
 , rho_spec(r), prod_spec(r);

rho = 1.4;
*rho = 1.7;

betalab(i,r) = sum(f$f_Lab(f),beta(f,i,r));

lab_sh("UnskLab",i,r)=fact_sh("UnskLab",i,r);
lab_sh("SkLab",i,r)=fact_sh("SkLab",i,r);

rho_spec(r) = rho - (rho-1) * sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * beta("SkLab",i,r) / betalab(i,r) );
prod_spec(r) = sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * beta("SkLab",i,r) / betalab(i,r) );

display rho_spec, prod_spec;

* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------
* DEFINE GENERAL EQUILIBRIUM MODEL
* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------

* the targetting is solved iteratively, in the following steps:
* 1) approximate by productivity
* 2) adjust trade / supplier costs
* 3) solve for trade costs
* 4) solve GE
* 5) re-adjust productivity, then reajust trade cost, solve GE, etc..

parameters      Sqr_stat_rep(TARG,T)
                Abs_stat_rep(TARG,T)
                chg_income_rep(TARG,T,r) per capita income
                chg_lambda_rep(TARG,T,r) lagrance multiplier
                chg_fd_rep(TARG,T,i,r) final demand
                chg_fd_qty_rep(TARG,T,i,r) final demand
                chg_absorption_rep(TARG,T,i,r) absorption
                chg_production_rep(TARG,T,i,r) production
                chg_production_qty_rep(TARG,T,i,R) production
                chg_Index_rep(TARG,T,i,r) Price index P
                chg_supplier_rep(TARG,T,i,r) Supplier effect S
                chg_fprice_rep(TARG,T,f,r) factor price
                chg_skillprem_rep(TARG,T,r) Skill premium
                chg_skillprem_last(r) Skill premium
                chg_income_ini(r) per capita income
                chg_fprice_ini(f,r) factor price
                chg_skillprem_ini(r) Skill premium;


* endogenous adjustment parameters for trade costs and productivity:
parameter prod_shock, bitc_shock;
prod_shock(r) = 1;
bitc_shock(r,s) = 1;


positive variables
                chg_income(r) per capita income
                chg_lambda(r) lagrance multiplier
                chg_fd(i,r) final demand
                chg_absorption(i,r) absorption
                chg_production(i,r) production
                chg_Index(i,r) Price index P
                chg_supplier(i,r) Supplier effect S
                chg_fprice(f,r) factor price
                chg_wage(f_Lab,r) wages
                ;


equations eq_demand, eq_budget, eq_absorption, eq_tradeprod, eq_supplier, eq_index, eq_wage;

* CLM preferences:
eq_demand(i,r)..     chg_fd(i,r) =e= chg_income(r)**sigma_p * chg_lambda(r)**(sigma_p-epsi(i))  * chg_index(i,r)**(1-sigma_p);

eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) *chg_fd(i,r));

eq_absorption(j,r).. chg_absorption(j,r) =e= gosh_f(j,r) *chg_fd(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r) *chg_production(i,r));

* Combining Trade and production:
eq_tradeprod(i,r)..  chg_production(i,r) =e=  sum(s,prod_sh(i,r,s) * chg_supplier(i,r) * bitc_shock(r,s)**(-theta(i)) *
                         chg_index(i,s)**theta(i) * chg_absorption(i,s));

* note: index = phi**(-1/theta) = P
eq_index(i,s)..      chg_index(i,s)**theta(i) * sum(r,trade_sh(i,r,s)*chg_supplier(i,r)*bitc_shock(r,s)**(-theta(i)) ) =e= 1;
* check for problems with zero-sum shares

eq_supplier(i,r)..   chg_supplier(i,r)  =e= prod(f$(f_noLab(f) and beta(f,i,r)),(chg_fprice(f,r)/prod_shock(r))**beta(f,i,r))**(-theta(i)) *
                         sum(f$(f_Lab(f) and beta(f,i,r)),(beta(f,i,r)/betalab(i,r))*(chg_fprice(f,r)/prod_shock(r))**(1-rho))**(-theta(i)*betalab(i,r)/(1-rho)) *
                         prod(j,chg_index(j,r)**intersh(j,i,r))**(-theta(i));

eq_wage(f_Lab,r)..   chg_wage(f_Lab,r)**rho =e= sum(i, lab_sh(f_Lab,i,r) * chg_production.L(i,r) /
                     (beta("UnskLab",i,r)/betalab(i,r)*chg_wage("UnskLab",r)**(1-rho) +  beta("SkLab",i,r)/betalab(i,r)*chg_wage("SkLab",r)**(1-rho)) ) ;


* define initial bounds:
chg_income.LO(r) =  0.0010;
chg_lambda.LO(r) =  0.001;
chg_fd.LO(i,r) =  0.001;
chg_absorption.LO(i,r) =  0.001;
chg_production.LO(i,r) =  0.001;
chg_Index.LO(i,r) =  0.001;
chg_supplier.LO(i,r) =  0.001;
chg_fprice.LO(f,r) =  0.001;

chg_income.UP(r) = 1000;
chg_lambda.UP(r) = 1000;
chg_fd.UP(i,r) = 1000;
chg_absorption.UP(i,r) = 1000;
chg_production.UP(i,r) = 1000;
chg_Index.UP(i,r) = 1000;
chg_supplier.UP(i,r) = 1000;
chg_fprice.UP(f,r) = 1000;

* define initial values:
chg_income.L(r) = 1;
chg_lambda.L(r) = 1;
chg_fd.L(i,r) = 1;
chg_absorption.L(i,r) = 1;
chg_production.L(i,r) = 1;
chg_Index.L(i,r) = 1;
chg_supplier.L(i,r) = 1;
chg_fprice.L(f,r) = 1;


chg_supplier.L(i,r) = prod(f$(f_noLab(f) and beta(f,i,r)),(chg_fprice.L(f,r)/prod_shock(r))**beta(f,i,r))**(-theta(i)) *
                         sum(f$(f_Lab(f) and beta(f,i,r)),(beta(f,i,r)/betalab(i,r))*(chg_fprice.L(f,r)/prod_shock(r))**(1-rho))**(-theta(i)*betalab(i,r)/(1-rho)) *
                         prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i));

* ------------------------
* real income calibration:
prod_shock(r) = income_target(r);

$if "%prodtarg%" == "load"  prod_shock(r) = reporting_("final","prod shock",r);
* small adjustment because ZMB sometimes causes model to crash.
$if "%prodtarg%" == "load"  prod_shock(r)$(sameas(r,"ZMB") and prod_shock(r) > 1.17)=  ((prod_shock(r)-1)*0.7) +1;

* ----------
* SUB-MODEL: initial solve for trade costs

variable tcost_(r);
positive variable  chg_btrade(i,r,s),chg_exports(r)  ;
equations eq_tcost_target, eq_exports, eq_btrade, eq_prod, eq_index2;


* targetting export/GDP
eq_tcost_target(r)..  tcost_target(r) * sum(i,production(i,r)/ sum((i.local),production(i,r)) * chg_production(i,r))
                        =e= chg_exports(r);

eq_exports(r).. chg_exports(r) =e= sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s)
                 / sum(i.local,sum(s.local$(not sameas(s,r)), trade("data",i,r,s))) *
                chg_btrade(i,r,s)));


eq_btrade(i,r,s)..  chg_btrade(i,r,s) =e= chg_supplier(i,r) * (1+ tcost_(r)$(not sameas(r,s)))**(-theta(i)) *
                         chg_index(i,s)**theta(i) * chg_absorption(i,s);

eq_prod(i,r)..  chg_production(i,r) =e=  sum(s,prod_sh(i,r,s) * chg_btrade(i,r,s));

eq_index2(i,s)..  chg_index(i,s)**theta(i) * sum(r,trade_sh(i,r,s)*chg_supplier(i,r)*
                                (1+ tcost_(r)$(not sameas(r,s)))**(-theta(i))) =e= 1;

chg_absorption.fx(i,r) = 1;
*chg_Index.fx(i,r) = 1;
chg_supplier.fx(i,r) = chg_supplier.L(i,r);

chg_exports.L(r) = 1;
chg_production.L(i,r) = 1;
chg_btrade.L(i,r,s) = 1;
tcost_.L(r) =0;

model tcost_targeting / eq_tcost_target.tcost_, eq_exports.chg_exports, eq_btrade.chg_btrade, eq_prod.chg_production, eq_index2.chg_index /;
tcost_targeting.iterlim = 1000000;

$if "%tradetarg%" == "yes"  solve tcost_targeting using MCP;
chg_absorption.LO(i,r) = 0;
chg_Index.lO(i,r) = 0;
chg_supplier.lO(i,r) = 0;
chg_absorption.UP(i,r) = 1e10;
chg_Index.UP(i,r) = 1e10;
chg_supplier.UP(i,r) = 1e10;

display tcost_.l, tcost_target;

* ------------------------
* feed these trade costs to the GE model and solve:

* transform into the trade shock entering GE model:
bitc_shock(r,s) = (1+ tcost_.L(r));
bitc_shock(r,r) = 1;

$if "%tradetarg%" == "load" bitc_shock(r,s) = reporting_("final","trade cost",r);


bitc_shock(r,r) = 1;

display bitc_shock;

* convergence parameter:
parameter        count, Sqr_stat, Abs_stat,
                target_conv target convergence param;
count = 0;
Abs_stat = 100;
Sqr_stat = 100;

* Reporting for initial iteration :
reporting("ini","export/GDP target",r) = tcost_target(r);

reporting("ini","trade cost",r) = (1+ tcost_.L(r));

reporting("ini","exports chg",r) = sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s) * ( chg_supplier.L(i,r) * (1+ tcost_.L(r)$(not sameas(r,s)))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))  ))
                         / sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s)));


reporting("ini","production chg",r) = sum(i, production(i,r)/ sum((i.local),production(i,r))
                        * sum(s,prod_sh(i,r,s) *
                        ( chg_supplier.L(i,r) * (1+ tcost_.L(r)$(not sameas(r,s)))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))));

reporting("ini","export/GDP chg",r) = reporting("ini","exports chg",r) / reporting("ini","production chg",r) ;


target_conv("ini","trade") = sum(r, reporting("ini","export/GDP chg",r) -  tcost_target(r));



* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------
* COMPUTE PARTIAL EQUILIBRIUM APPROXIMATIONS OF SKILL PREMIUM

* clm diffs
parameters incelast, epsibar_diff, chg_fd_approx, chg_absorption_approx, chg_fprice_justfd, chg_fprice_approx, chg_skillprem_approx, chg_skillprem_justfd;

*incelast(i,r) = sigma(i) * sum(j, fd(j,r)) / sum(j, sigma(j) * fd(j,r));
*chg_fd_approx(i,r) = 1 + (incelast(i,r) - 1) * (prod_shock(r) - 1)* target_rescale_factor;



epsibar_diff(r) = sum(j, epsilon_diff(j) * fd(j,r)) / sum(j, fd(j,r));

*incelast(i,r) = sigma(i) * sum(j, fd(j,r)) / sum(j, sigma(j) * fd(j,r));
incelast(i,r) = 1 + (1 - sigma_p) * (epsilon_diff(i) - epsibar_diff(r)) / epsibar_diff(r) ;
* CHECK THIS! (NOT USED FOR NOW)

display incelast, epsibar_diff;
* note: best to not rescale here (some negative values)
chg_fd_approx(i,r) = 1 + (incelast(i,r) - 1) * (prod_shock(r) - 1);

chg_fd_approx(i,r)$(chg_fd_approx(i,r) < 0.001) = 0.0001;


* Inverting Gosh matrix:

equations eq_absorption_auatarky;
* absorption = production:
eq_absorption_auatarky(j,r).. chg_absorption(j,r) =e= gosh_f(j,r) *chg_fd_approx(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r) * chg_absorption(i,r));

model Invert_gosh / eq_absorption_auatarky.chg_absorption /;
Invert_gosh.iterlim = 1000000;

solve Invert_gosh using MCP;

chg_absorption_approx(i,r) = chg_absorption.L(i,r);

chg_absorption_approx(i,r)$(not chg_absorption_approx(i,r)) = 1;


chg_fprice_justfd(f,r) = sum(i,fact_sh(f,i,r) * chg_fd_approx(i,r));
chg_fprice_approx(f,r) = sum(i,fact_sh(f,i,r) * chg_absorption_approx(i,r));

* SKILL PREMIUM APPROXIMATIONS USED IN THE PAPER:
chg_skillprem_approx(r) = sum(i,fact_sh("SkLab",i,r) * chg_absorption_approx(i,r)) / sum(i,fact_sh("UnskLab",i,r) * chg_absorption_approx(i,r));
chg_skillprem_justfd(r) = sum(i,fact_sh("SkLab",i,r) * chg_fd_approx(i,r)) / sum(i,fact_sh("UnskLab",i,r) * chg_fd_approx(i,r));


parameter decomposition;

decomposition(r,"Approx ratio") = sum(i,fact_sh("SkLab",i,r) * chg_absorption_approx(i,r)) / sum(i,fact_sh("UnskLab",i,r) * chg_absorption_approx(i,r));
decomposition(r,"JustFD approx ratio") = sum(i,fact_sh("SkLab",i,r) * chg_fd_approx(i,r)) / sum(i,fact_sh("UnskLab",i,r) * chg_fd_approx(i,r));
decomposition(r,"Approx diff") = sum(i,(fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * log(chg_absorption_approx(i,r)));
decomposition(r,"JustFD approx diff") = sum(i,(fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * log(chg_fd_approx(i,r)));
decomposition(r,"JustFD approx eps") = sum(i,(fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * incelast(i,r));
decomposition(r,"logpci") = log(pcexp(r));

parameter fact_sh_av;

fact_sh_av("SkLab",i) = sum(r,fact_sh("SkLab",i,r)) / sum(r,1);
fact_sh_av("UnskLab",i) = sum(r,fact_sh("UnskLab",i,r)) / sum(r,1);

decomposition(r,"Av shares - with IO") = sum(i,(fact_sh_av("SkLab",i)-fact_sh_av("UnskLab",i)) * log(chg_absorption_approx(i,r)));
decomposition(r,"Av shares - just FD") = sum(i,(fact_sh_av("SkLab",i)-fact_sh_av("UnskLab",i)) * incelast(i,r));

* also need to compute changes in absorption associated with changes in FD computed using average incelast
parameter chg_fd_approx_avgincel, chg_absorption_approx_avgincel;
chg_fd_approx_avgincel(i,r) = 1 + (   sum(s,incelast(i,s))/ sum(r.local,1)  - 1) * (prod_shock(r) - 1);

variable chg_absorption_avgincel;
equations eq_absorption_auatarky_avgincel;
* absorption = production:
eq_absorption_auatarky_avgincel(j,r).. chg_absorption_avgincel(j,r) =e= gosh_f(j,r) *chg_fd_approx_avgincel(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r) * chg_absorption_avgincel(i,r));

model Invert_gosh_avgincel / eq_absorption_auatarky_avgincel.chg_absorption_avgincel /;
Invert_gosh.iterlim = 1000000;

solve Invert_gosh_avgincel using MCP;

chg_absorption_approx_avgincel(i,r) = chg_absorption_avgincel.L(i,r);

chg_absorption_approx_avgincel(i,r)$(chg_absorption_approx_avgincel(i,r) < 0.001) = 0.001;

* USED IN PAPER:
decomposition(r,"Av incelast - with IO") = sum(i,(fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * log(chg_absorption_approx_avgincel(i,r))  );
decomposition(r,"Av incelast - just FD") = sum(i,(fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * sum(s,incelast(i,s))/ sum(r.local,1) );


* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------
* NOW SOLVE THE GENERAL EQUILIBRIUM MODEL

* FOUR STEPS PER LOOP:
* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
* 4) MANUALLY ADJUSTING FACTOR PRICES


* MODELS:

* 1) MODEL TO SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
model counterfactMCP_PE_justprices / eq_supplier.chg_supplier, eq_index.chg_index /;
counterfactMCP_PE_justprices.iterlim = 1000000;

* 2) MODEL SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
model counterfactMCP_PE_demand / eq_demand.chg_fd, eq_budget.chg_lambda /;
counterfactMCP_PE_demand.iterlim = 1000000;

* 3) MODEL FOR SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
model counterfactMCP_PE_fixedprice / eq_absorption.chg_absorption, eq_tradeprod.chg_production/;
counterfactMCP_PE_fixedprice.iterlim = 1000000;

* 4) MODEL FOR SOLVING FOR WAGES, TAKING PRODUCTION AS GIVEN:
model counterfactMCP_PE_wage / eq_wage.chg_wage/;
counterfactMCP_PE_wage.iterlim = 1000000;

***

* FULL GE for simultaneous solve:
*model counterfactMCP / eq_demand.chg_index, eq_budget.chg_lambda, eq_absorption.chg_absorption, eq_tradeprod.chg_production,
*eq_supplier.chg_supplier, eq_index.chg_fd, eq_fmc.chg_fprice, eq_income.chg_income /;


* INITIAL ADJUSTMENTS:
* Prices changes equal the income shock when we neglect changes in factor prices.
* just need to solve for change in demand and production: 3)

chg_fprice.FX(f,r) = chg_fprice_approx(f,r);
chg_Income.FX(r) = sum(f,endow_sh(f,r) * chg_fprice.L(f,r));

* In the income growth counterfactual, this is already the solution:
chg_supplier.L(i,r) = prod(f$beta(f,i,r),prod_shock(r)**beta(f,i,r))**theta(i) * prod(j,prod_shock(r)**intersh(j,i,r))**theta(i);

chg_index.L(i,r) = 1/prod_shock(r);

solve counterfactMCP_PE_justprices using MCP ;

chg_supplier.FX(i,r) = chg_supplier.L(i,r);
chg_index.FX(i,r) = chg_index.L(i,r);

solve counterfactMCP_PE_demand using MCP ;

chg_fd.FX(i,r) = chg_fd.L(i,r);
chg_lambda.FX(r) = chg_lambda.L(r);

solve counterfactMCP_PE_fixedprice using MCP ;

* Storing first iteration:
chg_fprice_ini(f,r) = sum(i,fact_sh(f,i,r) * chg_production.L(i,r));
chg_income_ini(r) = sum(f,endow_sh(f,r) * chg_fprice_ini(f,r));
chg_skillprem_ini(r) = sum(i,fact_sh("SkLab",i,r) * chg_production.L(i,r)) / sum(i,fact_sh("UnskLab",i,r) * chg_production.L(i,r));


* compute real income change and compare to target
reporting("initial","real income geo avg chg",r) = exp( sum((i), dem_sh(i,r)  * log(chg_fd.L(i,r)   /chg_Index.L(i,r)))  )     ;
reporting("initial","prod shock",r) = prod_shock(r);


display reporting;
parameter income_dev(r) income dev in pct terms;

* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------
* START TRADE/PRODUCTIVITY TARGETTING LOOPS HERE
loop(TARG,

* first : ------------
* re-adjust productivity shock


reporting(TARG,"real income geo avg chg",r) = exp( sum((i), dem_sh(i,r)  * log(chg_fd.L(i,r)   /
chg_Index.L(i,r)))  )     ;


* defining income dev as ratio. Probably best but troubled if close to 1.
income_dev(r)$(income_target(r) -1) = ( income_target(r)- reporting(TARG,"real income geo avg chg",r)) / (income_target(r) -1);

reporting(TARG,"income dev",r) = income_dev(r);

* then adjust prod_shock, using arbitrary adjustment rule:
* note: adjustment depends on whether prod_shock is currently above or below 1

* normal case:
$if "%prodtarg%" == "yes" prod_shock(r)$(prod_shock(r) >1) =1+ ((prod_shock(r)-1) * (1+ income_dev(r)/3));

* if goes below 1, flip target
$if "%prodtarg%" == "yes" prod_shock(r)$(prod_shock(r) le 1) =1+ ((prod_shock(r)-1) * (1 - income_dev(r)/3));

* put some bounds on adjustment (non-binding in equilibrium):
prod_shock(r)$(prod_shock(r) < 0.5) =0.5;
prod_shock(r)$(prod_shock(r) > 1.7) =1.7;

reporting(TARG,"prod shock",r) = prod_shock(r);

display income_dev, prod_shock;

* re-compute chg_supplier:
chg_supplier.fx(i,r) = prod(f$(f_noLab(f) and beta(f,i,r)),(chg_fprice.L(f,r)/prod_shock(r))**beta(f,i,r))**(-theta(i)) *
                         sum(f$(f_Lab(f) and beta(f,i,r)),(beta(f,i,r)/betalab(i,r))*(chg_fprice.L(f,r)/prod_shock(r))**(1-rho))**(-theta(i)*betalab(i,r)/(1-rho)) *
                         prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i));


***
* second : -------------
* re-adjust trade cost by re-solving Tcost model

* fixed inputs: chg_supplier, absorption, index
* endog: tcost_

chg_absorption.fx(i,r) = chg_absorption.l(i,r);


$if "%tradetarg%" == "yes" solve tcost_targeting using MCP;

chg_absorption.LO(i,r) = 0;
chg_Index.lO(i,r) = 0;
chg_supplier.lO(i,r) = 0;
chg_absorption.UP(i,r) = 1e10;
chg_Index.UP(i,r) = 1e10;
chg_supplier.UP(i,r) = 1e10;

* transform into the trade shock entering GE model:
bitc_shock(r,s) = (1+ tcost_.L(r));
$if "%tradetarg%" == "load" bitc_shock(r,s) = reporting_("final","trade cost",r);
bitc_shock(r,r) = 1;

* turn off LST file creastion (to avoid excessively large files)
OPTION SOLPRINT = OFF;

* ---------------------------------------------------------------------------
* ---------------------------------------------------------------------------
* START GENERAL EQUILIBRIUM LOOPS HERE

loop(T,

*for counterfactual in tau:
If ((Abs_stat > 0.0000001),

*while((Abs_stat > 0.0001 and count < 3),

* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:

chg_supplier.LO(i,r) = 0.001;
chg_supplier.UP(i,r) = 1000;
chg_index.LO(i,r) = 0.001;
chg_index.UP(i,r) = 1000;

solve counterfactMCP_PE_justprices using MCP ;

chg_supplier.FX(i,r) = chg_supplier.L(i,r);
chg_index.FX(i,r) = chg_index.L(i,r);


* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:

chg_fd.LO(i,r) = 0.001;
chg_fd.UP(i,r) = 1000;
chg_lambda.LO(r) = 0.001;
chg_lambda.UP(r) = 1000;

solve counterfactMCP_PE_demand using MCP ;

chg_fd.FX(i,r) = chg_fd.L(i,r);
chg_lambda.FX(r) = chg_lambda.L(r);


* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:

solve counterfactMCP_PE_fixedprice using MCP ;

* 4) SOLVING FOR WAGES, TAKING PRODUCTION AS GIVEN:

chg_wage.LO("SkLab",r) = 0.001;
chg_wage.UP("SkLab",r) = 1000;
chg_wage.LO("UnskLab",r) = 0.001;
chg_wage.UP("UnskLab",r) = 1000;

chg_production.FX(i,r) = chg_production.L(i,r);
chg_wage.L("SkLab",r) = chg_fprice.L("SkLab",r);
chg_wage.L("UnskLab",r) = chg_fprice.L("UnskLab",r);


solve counterfactMCP_PE_wage using MCP ;

chg_production.LO(i,r) = 0.001;
chg_production.UP(i,r) = 1000;

***

* 5) CONVERGENCE STATS + MANUALLY ADJUSTING FACTOR PRICES:

Sqr_stat = ( sum(r, sqr(1-chg_wage.L("UnskLab",r)/chg_fprice.L("UnskLab",r)) + sqr(1- chg_wage.L("SkLab",r)/chg_fprice.L("SkLab",r)) ) /  sum((r,f_Lab),1)  )**(1/2);
Abs_stat = sum(r, abs(1-chg_wage.L("UnskLab",r)/chg_fprice.L("UnskLab",r)) + abs(1-chg_wage.L("SkLab",r)/chg_fprice.L("SkLab",r)) ) /  sum((r,f_Lab),1);

count = count + 1;

display Sqr_stat, Abs_stat, count, chg_fprice.L;

* Reporting:
Sqr_stat_rep(TARG,T) = Sqr_stat;
Abs_stat_rep(TARG,T) =  Abs_stat;
chg_lambda_rep(TARG,T,r) = chg_lambda.L(r);
chg_fd_rep(TARG,T,i,r) = chg_fd.L(i,r);
chg_fd_qty_rep(TARG,T,i,r) = chg_fd.L(i,r) / chg_Index.L(i,r);
chg_absorption_rep(TARG,T,i,r) = chg_absorption.L(i,r);
chg_production_rep(TARG,T,i,r) = chg_production.L(i,r);
chg_production_qty_rep(TARG,T,i,R) = chg_production.L(i,r) / chg_Index.L(i,r);
chg_Index_rep(TARG,T,i,r) = chg_Index.L(i,r);
chg_supplier_rep(TARG,T,i,r) = chg_supplier.L(i,r);
chg_fprice_rep(TARG,T,f,r)$f_noLab(f) = sum(i,fact_sh(f,i,r) * chg_production.L(i,r));
chg_fprice_rep(TARG,T,"UnskLab",r) = chg_wage.L("UnskLab",r);
chg_fprice_rep(TARG,T,"SkLab",r) = chg_wage.L("SkLab",r);
chg_income_rep(TARG,T,r) = sum(f,endow_sh(f,r) * chg_fprice_rep(TARG,T,f,r));
chg_skillprem_rep(TARG,T,r) = chg_fprice_rep(TARG,T,"SkLab",r)/chg_fprice_rep(TARG,T,"UnskLab",r);
chg_skillprem_last(r) =       chg_fprice_rep(TARG,T,"SkLab",r)/chg_fprice_rep(TARG,T,"UnskLab",r);

* Adjusting by less for slow but smooth convergence for TRADE counterfactual:
chg_fprice.FX("SkLab",r)   = (4/5) * chg_fprice.L("SkLab",r)   + (1/5) * chg_wage.L("SkLab",r);
chg_fprice.FX("UnskLab",r) = (4/5) * chg_fprice.L("UnskLab",r) + (1/5) * chg_wage.L("UnskLab",r);
chg_fprice.FX("Land",r)    = (4/5) * chg_fprice.L("Land",r)    + (1/5) * sum(i,fact_sh("Land",i,r) * chg_production.L(i,r));
chg_fprice.FX("Capital",r) = (4/5) * chg_fprice.L("Capital",r) + (1/5) * sum(i,fact_sh("Capital",i,r) * chg_production.L(i,r));
chg_fprice.FX("NatlRes",r) = (9/10)* chg_fprice.L("NatlRes",r) + (1/10)* sum(i,fact_sh("NatlRes",i,r) * chg_production.L(i,r));

* Adjusting by less for slow but smooth convergence:
*chg_fprice.FX(f,r) = (3/4) * chg_fprice.L(f,r) + (1/4) * sum(i,fact_sh(f,i,r) * chg_production.L(i,r));

chg_income.FX(r) = sum(f,endow_sh(f,r) * chg_fprice.L(f,r));

* end of "GE" model loop:
);
display count;
);

* back to trade cost targetting:

reporting(TARG,"export/GDP target",r) = tcost_target(r);
reporting(TARG,"real income target",r) = income_target(r) ;
reporting(TARG,"real income geo avg chg",r) = exp( sum((i), dem_sh(i,r)  * log(chg_fd.L(i,r)   /
chg_Index.L(i,r)))  );


reporting(TARG,"trade cost",r) = (1+ tcost_.L(r));
reporting(TARG,"exports chg",r) = sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s) * ( chg_supplier.L(i,r) * (bitc_shock(r,s))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))  ))
                         / sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s)));

reporting(TARG,"production chg",r) = sum(i, production(i,r)/ sum((i.local),production(i,r))
                        * sum(s,prod_sh(i,r,s) *
                        ( chg_supplier.L(i,r) * (bitc_shock(r,s))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))));

reporting(TARG,"export/GDP chg",r) = reporting(TARG,"exports chg",r) / reporting(TARG,"production chg",r) ;


target_conv(TARG,"trade") = sum(r, reporting(TARG,"export/GDP chg",r) -  tcost_target(r));

display reporting, target_conv;

* end of  targetting  loop:
);

* FINAL REPORTING
reporting("final","pop",r)$pcexp(r) = pop(r);
reporting("final","logpci",r)$pcexp(r) = log(pcexp(r));
reporting("final","chg_skilprem_last RESCALED",r) = (chg_skillprem_last(r) -1) * target_rescale_factor +1;
reporting("final","chg_income nominal",r) = chg_income.l(r);
reporting("final","export/GDP target",r) = tcost_target(r);
reporting("final","trade cost",r) = (1+ tcost_.L(r));
reporting("final","prod shock",r) = prod_shock(r);
reporting("final","exports chg",r) = sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s) * ( chg_supplier.L(i,r) * (bitc_shock(r,s))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))  ))
                         / sum(i, sum(s$(not sameas(s,r)), trade("data",i,r,s)));

reporting("final","production chg",r) = sum(i, production(i,r)/ sum((i.local),production(i,r))
                        * sum(s,prod_sh(i,r,s) *
                        ( chg_supplier.L(i,r) * (bitc_shock(r,s))**(-theta(i)) *
                         chg_index.L(i,s)**theta(i) * chg_absorption.L(i,s))));

reporting("final","export/GDP chg",r) = reporting("final","exports chg",r) / reporting("final","production chg",r) ;


reporting("final","real income target",r)  = income_target(r) ;


reporting("final","real income geo avg chg",r) = exp( sum((i), dem_sh(i,r)  * log(chg_fd.L(i,r)   /
chg_Index.L(i,r)))  );

reporting("final","real income chg",r) = sum((i), dem_sh(i,r)  * chg_fd.L(i,r) )  / sum((i), dem_sh(i,r)  * chg_Index.L(i,r));
reporting("final","real income chg2",r) = sum((i), dem_sh(i,r)  * chg_fd.L(i,r)   /  chg_Index.L(i,r));

reporting("final","rhotilde",r) = rho - (rho-1) * sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r)) * beta("SkLab",i,r) / betalab(i,r) );


target_conv("final","trade") = sum(r, reporting("final","export/GDP chg",r) -  tcost_target(r));

parameter skillpremium;
skillpremium(r) = reporting("final","chg_skilprem_last RESCALED",r) -1 ;

* ---------------------------------------------------------
* ---------------------------------------------------------
* END MODEL
* ---------------------------------------------------------


* ---------------------------------------------------------
* TRADE COST COUNTERFACTUAL DECOMPOSITION

alias(r, rloop);

parameter log_chg_skillprem_approx;

* first approx, without trying to approximat the change in production (but still linear approx):
log_chg_skillprem_approx(T,"no approx",r) = log( sum(i, fact_sh("SkLab",i,r) * chg_production.L(i,r) ))
                                         - log(sum(i, fact_sh("UnskLab",i,r) * chg_production.L(i,r) ));

*1- drop the log TC change in front
*2- multiply each share by the relevant calibrated change in border effects (in the new direct effect -- no need to change other components of the decomposition)
*Note: in the very last term,  used the ratio "prod_sh_adj(i,r,r)/trade_sh(i,r,r)" instead of D/Y. Should be equivalent.

* First component ("direct effect"):
log_chg_skillprem_approx(T,"direct-refined",r) = - log(0.99) *  sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                               * theta(i) * ( sum(s, prod_sh(i,r,s)*(1-trade_sh(i,r,s)) ) - 2*prod_sh(i,r,r)*(1-trade_sh(i,r,r)) ) );

log_chg_skillprem_approx(T,"direct_Feb2019 alt",r) = -  sum(i$production(i,r), (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                 * theta(i) * ( sum(s$(not sameas(r,s)) ,prod_sh(i,r,s) * log(bitc_shock(r,s)) )
                - sum(s$(not sameas(r,s)), trade_sh(i,s,r) *  log(bitc_shock(s,r))) *  absorption(i,r)/ production(i,r)));

log_chg_skillprem_approx(T,"direct_Feb2019",r) = -  sum(i$trade_sh(i,r,r) , (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                * theta(i) * ( sum(s$(not sameas(r,s)) ,prod_sh(i,r,s) * log(bitc_shock(r,s)) )
- sum(s$(not sameas(r,s)), trade_sh(i,s,r) *  log(bitc_shock(s,r))) *   prod_sh(i,r,r) /  trade_sh(i,r,r) ));


* Combine these two effects to get effect of final demand (2nd component):
log_chg_skillprem_approx(T,"lambda",r) =  sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                                         * (-sigma(i)) * sum(s, gosh_f(i,s) * prod_sh(i,r,s) * log(chg_lambda.L(s)) )  );

log_chg_skillprem_approx(T,"P2",r) = sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                                         * (1-sigma(i)) * sum(s,prod_sh(i,r,s) * gosh_f(i,s) * log(chg_Index.L(i,s))) );


* third component (IO and GE feedback):
log_chg_skillprem_approx(T,"IO",r)= sum(i, (fact_sh("SkLab",i,r)-fact_sh("UnskLab",i,r))
                                         * sum(s,prod_sh(i,r,s) * sum(j,gosh_m(i,j,s)*log(chg_production.L(j,s))) ));


log_chg_skillprem_approx(T,"MR_Feb2019",rloop) = sum(i$trade_sh(i,rloop,rloop), (fact_sh("SkLab",i,rloop)-fact_sh("UnskLab",i,rloop)) * (
                                                 sum(s, prod_sh(i,rloop,s)*log(chg_supplier.L(i,rloop)*chg_Index.L(i,s)**theta(i)) )
                                                 - prod_sh(i,rloop,rloop)/trade_sh(i,rloop,rloop) *
                                                 sum(r, trade_sh(i,r,rloop)*log(chg_supplier.L(i,r)*chg_Index.L(i,rloop)**theta(i)) )
                                                 ) );


log_chg_skillprem_approx(T,"CES_scaling",r) = rho_spec(r);


log_chg_skillprem_approx(T,"Resum",r) =  ( log_chg_skillprem_approx(T,"direct_Feb2019",r) +
                                           log_chg_skillprem_approx(T,"MR_Feb2019",r) +
                                           log_chg_skillprem_approx(T,"IO",r) +
                                           log_chg_skillprem_approx(T,"P2",r) +
                                           log_chg_skillprem_approx(T,"lambda",r)  ) / log_chg_skillprem_approx(T,"CES_scaling",r);

parameter scaling;
scaling("rhotilde",r) = rho_spec(r);
scaling("rho",r) = rho;
scaling("skillprem",r) = skillpremium(r);


* ---------------------------------------------------------
* Export results

* also SP estimates to merge then load to stata
execute_unload 'results/SPestimates_tomerge_%demandest%_CLM/CF_prod_%prodtarg%_trade_%tradetarg%_%spec%_CLM.gdx', skillpremium, log_chg_skillprem_approx, scaling, decomposition;
*execute_unload 'SPestimates_tomerge_%demandest%_highersigma/unif_prod_%prodtarg%_trade_%tradetarg%_%spec%.gdx', skillpremium, log_chg_skillprem_approx;


* also unload full results:
execute_unload 'results/%demandest%_CLM/CF_Skill_prod_%prodtarg%_trade_%tradetarg%_%spec%_CLM.gdx',
                reporting,
                sigma,
                theta,
                Sqr_stat_rep,
                Abs_stat_rep,
                chg_income_rep,
                chg_lambda_rep,
                chg_fd_rep ,
                chg_absorption_rep ,
                chg_production_rep ,
                chg_Index_rep ,
                chg_supplier_rep ,
                chg_fprice_rep,
                chg_fd_qty_rep,
                chg_production_qty_rep,
                chg_skillprem_rep,
                chg_income_ini,
                chg_fprice_ini,
                chg_skillprem_ini,
                pcexp,
                incelast,
                chg_fd_approx,
                chg_absorption_approx,
                chg_fprice_justfd,
                chg_fprice_approx,
                chg_skillprem_approx,
                chg_skillprem_justfd,
                target_conv
                log_chg_skillprem_approx
                decomposition;


