cd C:\Projekte\homepage

********************************************************************************
clear all
matrix drop _all
clear mata
capture log close 
capture set matsize 8000 
capture set more off
program drop _all
sca drop _all
********************************************************************************

********************************************************************************
log using stata_cppml, replace
use stata_cppml_data, clear /* use data from a Monte Carlo run */
********************************************************************************

********************************************************************************
*** Variables
********************************************************************************
label var ex "exporter"
label var im "importer"
label var y  "observed trade flow"
label var yi "gross production share" 
label var yj "expediture share"
label var V  "missingess indicator"
label var x1 "border dummy 1 if ex ~=im"
label var x2 "log distance"
label var V  "missingness indicator"
********************************************************************************

********************************************************************************
*** Counterfactual *** 
********************************************************************************
gen x1c= 0 /* no borders */
label var x1c "counterfactual no border"

********************************************************************************
*** Globals *** 
********************************************************************************
global K=2
global b=20
global groups=2 /* country groups for counterfactuals */
global sig=5 /*elasticity of substitution*/
global re =  "x1  x2" 
global recf= "x1c x2" 
global dum=  "d1 d2 d3 d4 d5 d6" /* Dummies for breakdown of trade flows*/
global row_dmpcf="dom-small dom-large small-small large-large small-large large-small"
global row_dw="dom-small dom-large"
********************************************************************************

********************************************************************************
*** Unconstrained PPML ***
********************************************************************************
***ppmlhdfe y $re if V==1 , absorb(ex im) d	 nocons
glm y $re ibn.im ib$b.ex if V==1, nocons   ///
	 irls  robust  family(poisson)  // Plain unconstraine PPML  ib$b.ex   ibn.im
predict ypu 

/*
*** check addendum (i)
gen yppu=y if V==1
replace yppu=ypu if V==0
glm yppu $re ibn.im ib$b.ex, nocons   ///
	 irls  robust  family(poisson)  // Plain unconstraine PPML  ib$b.ex   ibn.im
*/
 
mat b0=e(b)
mat bppml=b0'
mat Vppml=e(V)
mat Vppml=Vppml[1..2,1..2]
mat Sppml=diag(vecdiag(Vppml))
mat Sppml=cholesky(Sppml)
mat bppml=bppml[1..2,.]
mat tppml=inv(Sppml)*bppml
mat Sppml=vecdiag(Sppml)'
mat out_ppml=(bppml, Sppml, tppml)
********************************************************************************

********************************************************************************
*** Initialize loop *** 
********************************************************************************
sca tol=0.000000000001 /*tolerance for convergence*/
sca itol=10000
sca iter=1

gen double thxm=yi*yj /*dependent variable for solver */
gen double sf=0
gen double za=0 
gen double r=0

qui foreach re of global re {
replace  za=za+_b[`re']*`re'
}
********************************************************************************

********************************************************************************
*** Loop starting for CPPML estimator *** 
********************************************************************************
qui while itol > tol & iter < 200 {
capture drop m

*** Inner loop starting solve for equilibrium *** 
*** glm thxm  ibn.im ib$b.ex , nocons offset(za) irls   family(poisson)
ppmlhdfe thxm , absorb(ex im) offset(za) d	 nocons

predict double m
replace  r=(y-m)/m
replace sf=za+(y-m)/m
replace sf=za if V==0

*** Iteration step in outer loop *** 
*** qui reg sf $re ibn.im ib$b.ex [aweight=m], nocons robust 
reghdfe sf $re  [aweight=m],  vce(robust) absorb(ex im) nocons

replace za=0,  nopromote 
foreach re of global re {
replace za=za+_b[`re']*`re'   /*update trade costs za*/
}

*** Prepare for next iteration step ***
qui {
sca iter=iter+1
mat b=e(b)
mat diff= mreldif(b,b0)
sca itol=diff[1,1]
mat b0=b /*update b */
}
if 10*int(iter/10)==iter {
di as red iter " "  itol
}
}
********************************************************************************

********************************************************************************
*** z'alpha, base and counterfactual
********************************************************************************
gen zacf=0
gen temp=x1
drop x1 
ren x1c x1 
foreach re of global re  {
replace    zacf=zacf+_b[`re']*`re'   
}
ren x1 x1c
ren temp x1
ren za zaba

label var m "prediced trade flow"
label var x1 "border"
label var x1c "counterfactual no border"
********************************************************************************

********************************************************************************
*** Calculate standard errors of CPPML estimates
********************************************************************************
qui  mata 
bcppml=st_matrix("b")'
y=st_data(., "y") 
m=st_data(., "m")
V=st_data(., "V")
Z=st_data(., "$re")
Dm=st_data(., "ibn.im")
Dx=st_data(., "ibn.ex")
Dx=Dx[., 1..$b-1]
D=(Dm,Dx)
M=diag(m)
V=diag(V)

Q=I(rows(Z))-M*D*luinv(D'*M*D)*D'
nob=J($b^2,1,1)'*V*J($b^2,1,1)
GIZZ=invsym(Z'*Q*V*M*V*Q'*Z)
Vcppml=(nob/(nob-1))*GIZZ*Z'*Q*V*(diag((y-m):*(y-m)))*V*Q'*Z*GIZZ'

Scppml=diagonal(Vcppml):^0.5
tcppml =bcppml[1..$K,.]:/ Scppml[1..$K,.]
out_cppml=(bcppml[1..$K,.], Scppml[1..$K,.], tcppml[1..$K,.] )
st_matrix("out_cppml", out_cppml)
end
********************************************************************************

********************************************************************************
*** Country-pair groups, large is fourth quartile
********************************************************************************
qui {
xtile sizex=yi, nq(4)
xtile sizem=yj, nq(4)

replace sizex=1 if sizex<4
replace sizex=2 if sizex==4
replace sizem=1 if sizem<4
replace sizem=2 if sizem==4

gen co = 1 if  sizem==1 & ex==im 
replace co = 2  if sizem==2 & ex==im 
replace co = 3 if  sizex==1 & sizem==1 & ex~=im
replace co = 4 if  sizex==2 & sizem==2 & ex~=im
replace co = 5 if  sizex==1 & sizem==2 & ex~=im
replace co = 6 if  sizex==2 & sizem==1 & ex~=im

label define lco 1 "dom-small", add
label define lco 2 "dom-large", add
label define lco 3 "small-small", add
label define lco 4 "large-large", add
label define lco 5 "small-large", add
label define lco 6 "large-small", add

label values co lco
label var co "country-pair groups"
}
********************************************************************************

********************************************************************************
*** Dummies for groups of bilateral flows to form S 
********************************************************************************
qui tab co , gen(d)	
qui forvalues i = 1(1)6 {
qui sum  d`i' 
qui replace d`i' = d`i'/r(sum)
}
********************************************************************************

********************************************************************************
*** Solve for counterfactual
********************************************************************************
qui glm thxm  ibn.im ib$b.ex , nocons offset(zacf) irls   family(poisson)
predict double phicf, xb
replace phicf= phicf-zacf
********************************************************************************

********************************************************************************
*** Delta method for counterfactual changes
********************************************************************************
qui mata 
printf("begin vcppmlcf \n")
Zcf=st_data(., "$recf")
S=st_data(., "$dum")
S=S'
zaba=st_data(., "zaba")
zacf=st_data(., "zacf")
phicf=st_data(., "phicf") 
mcf=exp(zacf+phicf)

M=diag(m)
Mcf=diag(mcf)
dmpcf =S*diagonal(  Mcf*luinv(M)-I($b^2) )
printf("end dmpcf  \n")

MI=luinv(M) 
MIcf=luinv(Mcf) 
Gacf=(Mcf-Mcf*D*luinv(D'*Mcf*D)*D'*Mcf)*Zcf
Ga=(M-M*D*(luinv(D'*M*D))*D'*M)*Z
Gacf=MIcf*Gacf
Ga=MI*Ga
printf("end some matrices \n")

Vdmpcf=S*MI*Mcf*(Gacf-Ga)*Vcppml*(Gacf-Ga)'*Mcf*MI*S'
Sdmpcf=diagonal(Vdmpcf:^0.5)
tdmpcf=dmpcf:/Sdmpcf
printf("done Vdmpcf \n")

dw=diagonal(MI*Mcf)
ddw=diag(  (1/(1-$sig)) *( dw :^ (1/(1-$sig)) ) )
dw=dw:^(1/(1-$sig))
printf("done dw ddw  \n")
Vdw=S*ddw*(Gacf-Ga)*Vcppml*(Gacf-Ga)'*ddw*S'
printf("done Vdw  \n")

dw=S*dw
dw=dw[1..$groups,.]
Vdw=Vdw[1..$groups,1..$groups]
Sdw=diagonal(Vdw:^0.5)
dw=dw- J(rows(dw),1,1)
tdw=dw:/Sdw

st_matrix("dw", dw)
st_matrix("Sdw", Sdw)
st_matrix("tdw", tdw)
st_matrix("dmpcf", dmpcf)
st_matrix("Vdmpcf", Vdmpcf)
st_matrix("Sdmpcf", Sdmpcf)
st_matrix("tdmpcf", tdmpcf)

printf("done vcppmlcf  \n")
end
drop zaba zacf sf 
********************************************************************************

********************************************************************************
*** Collect results in matrices
********************************************************************************
qui {
mat dmpcf=100*dmpcf
mat dmpcf_u= dmpcf-100*1.96*Sdmpcf
mat dmpcf_o= dmpcf+100*1.96*Sdmpcf
mat dw= 100*dw
mat dw_u= dw-100*1.96*Sdw
mat dw_o= dw+100*1.96*Sdw

mat out_dmp=(dmpcf, tdmpcf, dmpcf_u, dmpcf_o)
mat out_dw=(dw, tdw, dw_u , dw_o)

mat colnames out_dmp = dmp  t dmp_u dmp_o
mat colnames out_dw= dw t dw_u dw_o
mat rownames out_dmp = $row_dmpcf
mat rownames out_dw = $row_dw
mat rownames out_ppml = $re
mat colnames out_ppml = b s t
mat rownames out_cppml = $re
mat colnames out_cppml = b s t
}
********************************************************************************

********************************************************************************
*** Display results
********************************************************************************

*** Observed (y) and predicted (m) trade flows by country-pair group ***
label var ypu "predicted unconstrained PPML"
label var m   "predicted constrained PPML"

table ex, c(sum V sum y sum ypu sum m m yi) row format(%6.3f)   stubwidth(15) 
table im, c(sum V sum y sum ypu sum m m yj) row format(%6.3f)   stubwidth(15) 
table co, c(sum y sum ypu sum m) row format(%6.3f)   stubwidth(15)     

*** Estimation results ***
estout matrix(out_ppml,  fmt(3)),title("Unconstrained PPML estimation results")
estout matrix(out_cppml,  fmt(3)),title("CPPML estimation results")
estout matrix(out_dmp, fmt(2)), title("Counterfactual effects on trade flows in percent")
estout matrix(out_dw,  fmt(2)), title("Counterfactual welfare effects in percent")