Prof. Francis Petterini

Paper: https://link.springer.com/article/10.1007/s00181-023-02488-6

Data: http://petterini.paginas.ufsc.br/files/2023/06/data_private_policing.txt

Dofile:

/*
cls
clear all
set obs 10000

*set-up a grid:
matrix g = J(sqrt(_N),1,0)
local b = 3
local t = sqrt(_N)
local i = 1
while `i' <= `t' {
 matrix g[`i',1] = -`b'+(`i'-1)/((`t'-1)/(2*`b'))
local i = `i' + 1
}
matrix g1 = g#J(`t',1,1)
matrix g2 = J(`t',1,1)#g 
svmat double g1, name(x)
svmat double g2, name(y)
rename x1 v1
rename y1 v2

*set-up margins mean 0 sd 1:
gen G1_gauss = normal(v1)
gen g1_gauss = normalden(v1)
gen G2_gauss = normal(v2)
gen g2_gauss = normalden(v2)
scalar b = sqrt(6)/_pi
scalar a = digamma(1)*b
gen vv1 = (v1-a)/b
gen G1_gumbel = exp(-exp(-vv1))
gen g1_gumbel = exp(-vv1-exp(-vv1))/b
gen vv2 = (v2-a)/b
gen G2_gumbel = exp(-exp(-vv2))
gen g2_gumbel = exp(-vv2-exp(-vv2))/b
replace vv1 = (v1+a)/b
replace vv2 = (v2+a)/b
gen G1_r_gumbel = 1-exp(-exp(vv1))
gen g1_r_gumbel = exp(vv1-exp(vv1))/b
gen G2_r_gumbel = 1-exp(-exp(vv2))
gen g2_r_gumbel = exp(vv2-exp(vv2))/b

*comparing margins:
tw (line g1_gumbel v1, sort lcolor(red) lpattern(solid)) ///
 (line g1_gauss v1, sort lcolor(black) lpattern(solid)) ///
 (line g1_r_gumbel v1, sort lcolor(blue) lpattern(solid)), ///
 ytitle(density) xtitle("error with mean = 0 and sd = 1") xla(-3(1)3, axis(1)) xsize(28) ysize(10) scheme(s1mono) ///
 legend(cols(3) label(1 "Gumbel") label(2 "Gaussian") label(3 "reverse-Gumbel"))
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\marginals.eps", as(eps) replace

*gaussian with gaussian margins:
scalar theta = .707 /* tau = .5 */
gen C12 = (exp(-(invnormal(G1_gauss)^2-2*theta*invnormal(G1_gauss)*invnormal(G2_gauss)+invnormal(G2_gauss)^2)/(2*(1-theta^2)))/(2*_pi*sqrt(1-theta^2)))/ ///
 (normalden(invnormal(G1_gauss))*normalden(invnormal(G2_gauss)))
gen f12 = g1_gauss*g2_gauss*C12
replace f12 = 0 if f12 == .
label variable f12 "density"
tw (contour f12 v1 v2, ccolors(gs14) levels(20) zlabel(#5, format(%4.2f))), xtitle(v{subscript:1}) ytitle(v{subscript:2}) ///
 xla(-3(1)3, axis(1)) yla(-3(1)3, axis(1)) ysize(14) xsize(20) scheme(s1mono)
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\gaussian_gauss_gauss.png", as(png) replace

*clayton I:
scalar theta = 2 /* tau = .5 */
gen C = (G1_gauss^(-theta)+G2_gauss^(-theta)-1)^(-1/theta)
replace C12 = ((1+theta)*(C^(1+2*theta)))/((G1_gauss*G2_gauss)^(1+theta))
replace f12 = g1_gauss*g2_gauss*C12
replace f12 = 0 if f12 == .
label variable f12 "density"
tw (contour f12 v1 v2, ccolors(gs14) levels(20) zlabel(#5, format(%4.2f))), xtitle(v{subscript:1}) ytitle(v{subscript:2}) ///
 xla(-3(1)3, axis(1)) yla(-3(1)3, axis(1)) ysize(14) xsize(20) scheme(s1mono)
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\claytonI_gauss_gauss.png", as(png) replace

*clayton II:
gen D = ((1-G1_gauss)^(-theta)+(1-G2_gauss)^(-theta)-1)^(-1/theta)
replace C12 = ((1+theta)/(((1-G1_gauss)*(1-G2_gauss))^(1+theta)))*(D^(1+2*theta))
replace f12 = g1_gauss*g2_gauss*C12
replace f12 = 0 if f12 == .
label variable f12 "density"
tw (contour f12 v1 v2, ccolors(gs14) levels(20) zlabel(#5, format(%4.2f))), xtitle(v{subscript:1}) ytitle(v{subscript:2}) ///
 xla(-3(1)3, axis(1)) yla(-3(1)3, axis(1)) ysize(14) xsize(20) scheme(s1mono)
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\claytonII_gauss_gauss.png", as(png) replace

*gaussian with gumbel and reverser gumbel margins:
scalar theta = .707 /* tau = .5 */
replace C12 = (exp(-(invnormal(G1_gumbel)^2-2*theta*invnormal(G1_gumbel)*invnormal(G2_r_gumbel)+invnormal(G2_r_gumbel)^2)/(2*(1-theta^2)))/(2*_pi*sqrt(1-theta^2)))/ ///
 (normalden(invnormal(G1_gumbel))*normalden(invnormal(G2_r_gumbel)))
replace f12 = g1_gumbel*g2_r_gumbel*C12
replace f12 = 0 if f12 == .
label variable f12 "density"
tw (contour f12 v1 v2, ccolors(gs14) levels(20) zlabel(#5, format(%4.2f))), xtitle(v{subscript:1}) ytitle(v{subscript:2}) ///
 xla(-3(1)3, axis(1)) yla(-3(1)3, axis(1)) ysize(14) xsize(20) scheme(s1mono)
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\gaussian_gumbel_r_gumbel.png", as(png) replace
*/

cls
clear all
import delimited "E:\OneDrive - UFSC\artigo_policing\data\data.txt"

label variable id "city id"
label variable pop "number of inhabitants"
label variable police "police patrol officers"
label variable hom "firearm homicides"
label variable psec "private security guards"
label variable fsuicides "firearm suicides"
label variable suicides "all suicides"
label variable prison "prison industrial complex = 1"
label variable gdp "gdp per capita R$"
label variable gini "gini index 0-100"
label variable jobs "formal jobs %"
label variable unemp "unemployment %"
label variable poprural "rural population %"
label variable popmen "young men population %"

gen y1 = 1000*hom/pop
gen y2 = psec/police

*table 3:
gen y1_positive = y1 > 0
gen y2_positive = y2 > 0
tab y1_positive y2_positive

su gdp, meanonly 
replace gdp = (gdp-r(min))/(r(max)-r(min))
gen guns = cond(suicides==0,0,fsuicides/suicides)

foreach var of varlist gini jobs unemp poprural popmen {
 replace `var' = `var'/100
}
 
gen z2 = prison

/*
ktau y1 y2 z1 z2, stats(taub p)
*/

global X gini gdp jobs unemp poprural popmen guns

*table 4:
tabstat y1 y2 $X z2, statistics(mean median sd min max)

ivtobit y1 $X (y2 = z2) if y2 > 0, ll(0) first nolog
rivtest, ci points(500) gridmult(14)

ivregress 2sls y1 $X (y2 = z2)

ivregress 2sls y1 $X (y2 = z2) if y1 > 0 & y2 > 0

global X gini gdp jobs unemp poprural popmen

gen z1 = guns

tobit y1 $X z1 z2, ll(0)
matrix start1 = (e(b)[1..1,1..6],e(b)[1,9],e(b)[1,7],e(b)[1,8])
scalar sd1 = ln(sqrt(e(b)[1,10]))
tobit y2 $X z1 z2, ll(0)
matrix start2 = (e(b)[1..1,1..6],e(b)[1,9],e(b)[1,7],e(b)[1,8])
scalar sd2 = ln(sqrt(e(b)[1,10]))

*define start values:
matrix start = (start1,sd1,start2,sd2,0)
matrix coleq start = eq1 eq1 eq1 eq1 eq1 eq1 eq1 eq2 eq3 eq4 eq5 eq5 eq5 eq5 eq5 eq5 eq5 eq6 eq7 eq8 eq9
matrix colnames start = gdp gini jobs unemp poprural popmen _cons _cons _cons _cons gdp gini jobs unemp poprural popmen _cons _cons _cons _cons _cons
matrix list start

*define copula:
program define copula
 version 17.0
 args lf xb1 b11 b12 s1 xb2 b21 b22 s2 r 
 tempvar sigma1 sigma2 theta beta11 beta12 beta21 beta22 v1 v2 G1 G2 g1 g2 D C C1 C2 C12 f1P f2P f12
 quietly {
 local y1 "$ML_y1"
 local y2 "$ML_y2"
 local z1 "$ML_y3"
 local z2 "$ML_y4"
 local i "$ML_y5"
 local j "$ML_y6"
 local k "$ML_y7"
 *sigmas and theta definitions:
 gen double `sigma1' = exp(`s1') 
 gen double `sigma2' = exp(`s2')
 gen double `theta' = (2/_pi)*atan(`r')
 replace `theta' = exp(`r') if `k' != 1
 *gammas definitions:
 gen double `beta11' = `b11'
 gen double `beta12' = `b12'
 gen double `beta21' = `b21'
 gen double `beta22' = `b22' 
 *erros definitions:
 gen double `v1' = (`y1'-`beta11'*`z1'-`beta12'*`z2'-`xb1')/`sigma1'
 gen double `v2' = (`y2'-`beta21'*`z1'-`beta22'*`z2'-`xb2')/`sigma2'
 *marg1 definitions:
 gen double `G1' = normal(`v1')
 gen double `g1' = normalden(`v1')
 replace `G1' = exp(-exp(-`v1')) if `i' == 2
 replace `g1' = exp(-`v1'-exp(-`v1')) if `i' == 2
 replace `G1' = 1-exp(-exp(`v1')) if `i' == 3
 replace `g1' = exp(`v1'-exp(`v1')) if `i' == 3
 *marg2 definitions:
 gen double `G2' = normal(`v2')
 gen double `g2' = normalden(`v2') 
 replace `G2' = exp(-exp(-`v2')) if `j' == 2
 replace `g2' = exp(-`v2'-exp(-`v2')) if `j' == 2 
 replace `G2' = 1-exp(-exp(`v2')) if `j' == 3
 replace `g2' = exp(`v2'-exp(`v2')) if `j' == 3 
 *copula and its derivates:
 gen double `D' = binormal(invnormal(`G1'),invnormal(`G2'),`theta')
 replace `D' = (`G1'^(-`theta')+`G2'^(-`theta')-1)^(-1/`theta') if `k' == 2 /* 0 */
 *replace `D' = ((1-`G1')^(-`theta')+`G2'^(-`theta')-1)^(-1/`theta') if `k' == 2 /* 90 */
 replace `D' = ((1-`G1')^(-`theta')+(1-`G2')^(-`theta')-1)^(-1/`theta') if `k' == 3 /* 180 */
 *replace `D' = (`G1'^(-`theta')+(1-`G2')^(-`theta')-1)^(-1/`theta') if `k' == 3 /* 270 */
 gen double `C' = `D'
 replace `C' = `D' if `k' == 2 /* 0 */
 *replace `C' = `G2'-`D' if `k' == 2 /* 90 */
 replace `C' = `G1'+`G2'-1+`D' if `k' == 3 /* 180 */
 *replace `C' = `G1'-`D' if `k' == 3 /* 270 */
 gen double `C1' = normal((invnormal(`G2')-`theta'*invnormal(`G1'))/sqrt(1-`theta'^2)) 
 replace `C1' = (`D'/`G1')^(1+`theta') if `k' == 2 /* 0 */ 
 *replace `C1' = (`D'/(1-`G1'))^(1+`theta') if `k' == 2 /* 90 */ 
 replace `C1' = `g1'-(`D'/(1-`G1'))^(1+`theta') if `k' == 3 /* 180 */
 *replace `C1' = `g1'-(`D'/`G1')^(1+`theta') if `k' == 3 /* 270 */
 gen double `C2' = normal((invnormal(`G1')-`theta'*invnormal(`G2'))/sqrt(1-`theta'^2)) 
 replace `C2' = (`D'/`G2')^(1+`theta') if `k' == 2 /* 0 */
 *replace `C2' = `g2'-(`D'/`G2')^(1+`theta') if `k' == 2 /* 90 */
 replace `C2' = `g2'-(`D'/(1-`G2'))^(1+`theta') if `k' == 3 /* 180 */
 *replace `C2' = (`D'/(1-`G2'))^(1+`theta') if `k' == 3 /* 270 */ 
 gen double `C12' = (exp(-(invnormal(`G1')^2-2*`theta'*invnormal(`G1')*invnormal(`G2')+invnormal(`G2')^2)/(2*(1-`theta'^2)))/(2*_pi*sqrt(1-`theta'^2)))/ ///
 (normalden(invnormal(`G1'))*normalden(invnormal(`G2')))
 replace `C12' = ((1+`theta')/((`G1'*`G2')^(1+`theta')))*(`D'^(1+2*`theta')) if `k' == 2 /* 0 */
 *replace `C12' = ((1+`theta')/(((1-`G1')*`G2')^(1+`theta')))*(`D'^(1+2*`theta')) if `k' == 2 /* 90 */
 replace `C12' = ((1+`theta')/(((1-`G1')*(1-`G2'))^(1+`theta')))*(`D'^(1+2*`theta')) if `k' == 3 /* 180 */ 
 *replace `C12' = ((1+`theta')/((`G1'*(1-`G2'))^(1+`theta')))*(`D'^(1+2*`theta')) if `k' == 3 /* 270 */ 
 *terms of the likelihhod:
 gen double `f1P' = (`g1'/`sigma1')*`C1'
 gen double `f2P' = (`g2'/`sigma2')*`C2'
 gen double `f12' = (`g1'/`sigma1')*(`g2'/`sigma2')*`C12'
 *setup lf for reasonable bounds
 tempname minv
 scalar `minv' = 1e-20 
 replace `lf' = ln(min(1-`minv',max(`minv',`C'))) if `y1' == 0 & `y2' == 0
 replace `lf' = ln(max(`minv',`f1P')) if `y1' > 0 & `y2' == 0
 replace `lf' = ln(max(`minv',`f2P')) if `y1' == 0 & `y2' > 0
 replace `lf' = ln(max(`minv',`f12')) if `y1' > 0 & `y2' > 0
 }
end

*i and j define marginals, and k defines copula:
gen i = 1
gen j = 2
gen k = 2

/*

*checking:
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () ()
ml check
replace k = 2
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () ()
ml check
replace k = 3
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () ()
ml check
replace k = 1

*/

constraint 1 _b[eq9:_cons] = 0

*checking (2 tobit = these results):
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), constraint(1) technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml init start
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip

ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml init start
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
matrix start1 = e(b)

replace k = 2
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), constraint(1) technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml init start
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
matrix start2 = e(b)
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml init start2
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
matrix start2 = e(b)

/*

replace k = 3
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), constraint(1) technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
matrix start3 = e(b)
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), technique(nr 25 bhhh 5 dfp 5 bfgs 5 nr 10)
ml init start3
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
matrix start3 = e(b)

*procedure:
tic
local kk = 1
 while `kk' <= 3 {
 replace k = `kk'
 matrix L = J(3,3,0)
 matrix A = J(3,3,0) 
 matrix B11 = J(3,3,0) 
 matrix B12 = J(3,3,0) 
 matrix B21 = J(3,3,0) 
 matrix B22 = J(3,3,0) 
 matrix G1 = J(3,3,0) 
 matrix G2 = J(3,3,0) 
 matrix GG = J(3,3,0)
 matrix T = J(3,3,0)
 local ii = 1
 while `ii' <= 3 {
 replace i = `ii' 
 local jj = 1 
 while `jj' <= 3 {
 replace j = `jj'
 ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), technique(nr 20 bhhh 10 dfp 5 bfgs 5 nr 5 bhhh 5)
 ml init start`kk'
 ml max, difficult iterate(50) tolerance(1e-3) ltolerance(1e-3) nrtolerance(1e-3) skip 
 estimates store copula_`ii'`jj'_`kk'
 matrix L[`ii',`jj'] = e(ll)
 matrix A[`ii',`jj'] = 2*(e(k)-e(ll)) 
 matrix B11[`ii',`jj'] = _b[eq2:_cons]
 matrix B12[`ii',`jj'] = _b[eq3:_cons]
 matrix B21[`ii',`jj'] = _b[eq6:_cons]
 matrix B22[`ii',`jj'] = _b[eq7:_cons]
 matrix G1[`ii',`jj'] = B12[`ii',`jj']/B22[`ii',`jj']
 matrix G2[`ii',`jj'] = B21[`ii',`jj']/B11[`ii',`jj']
 matrix GG[`ii',`jj'] = G1[`ii',`jj']*G2[`ii',`jj']
 matrix T[`ii',`jj'] = _b[eq9:_cons]
 local jj = `jj' + 1
}
 local ii = `ii' + 1
}
 matrix L = vec(L')
 matrix A = vec(A')
 matrix B11 = vec(B11')
 matrix B12 = vec(B12')
 matrix B21 = vec(B21')
 matrix B22 = vec(B22')
 matrix G1 = vec(G1')
 matrix G2 = vec(G2')
 matrix GG = vec(GG')
 matrix T = vec(T')
 matrix copula_`kk' = (L,A,B11,B12,B21,B22,G1,G2,GG,T)
local kk = `kk' + 1
}
toc

*main results:
matrix list copula_1
matrix list copula_2
matrix list copula_3

*/

*reduced estimates:
replace i = 1
replace j = 2
replace k = 2
ml model lf copula (y1 y2 z1 z2 i j k = $X) () () () ($X) () () () (), technique(nr 20 bhhh 10 dfp 5 bfgs 5 nr 5 bhhh 5)
ml init start2
ml max, difficult iterate(50) tolerance(1e-5) ltolerance(1e-5) nrtolerance(1e-5) skip 
estimates store copula_best
matrix reduced = e(b)

*structural estimates:
scalar gamma1 = _b[eq3:_cons]/_b[eq7:_cons]
scalar gamma2 = _b[eq6:_cons]/_b[eq2:_cons]
scalar a11 = (1-gamma1*gamma2)*_b[eq2:_cons]
scalar a22 = (1-gamma1*gamma2)*_b[eq7:_cons]
scalar a1_cons = _b[eq1:_cons]-gamma1*_b[eq5:_cons]
scalar a1gdp = _b[eq1:gdp]-gamma1*_b[eq5:gdp]
scalar a1gini = _b[eq1:gini]-gamma1*_b[eq5:gini]
scalar a1jobs = _b[eq1:jobs]-gamma1*_b[eq5:jobs]
scalar a1unemp = _b[eq1:unemp]-gamma1*_b[eq5:unemp]
scalar a1poprural = _b[eq1:poprural]-gamma1*_b[eq5:poprural]
scalar a1popmen = _b[eq1:popmen]-gamma1*_b[eq5:popmen]
scalar a2_cons = _b[eq5:_cons]-gamma2*_b[eq1:_cons]
scalar a2gdp = _b[eq5:gdp]-gamma2*_b[eq1:gdp]
scalar a2gini = _b[eq5:gini]-gamma2*_b[eq1:gini]
scalar a2jobs = _b[eq5:jobs]-gamma2*_b[eq1:jobs]
scalar a2unemp = _b[eq5:unemp]-gamma2*_b[eq1:unemp]
scalar a2poprural = _b[eq5:poprural]-gamma2*_b[eq1:poprural]
scalar a2popmen = _b[eq5:popmen]-gamma2*_b[eq1:popmen]
scalar sigma1 = exp(_b[eq4:_cons])
scalar sigma2 = exp(_b[eq8:_cons])*(sqrt(6)/_pi)
scalar tau = exp(_b[eq9:_cons])/(exp(_b[eq9:_cons])+2)

*delta method:
nlcom (gamma1: _b[eq3:_cons]/_b[eq7:_cons]) (gamma2: _b[eq6:_cons]/_b[eq2:_cons]) ///
 (a11: (1-gamma1*gamma2)*_b[eq2:_cons]) (a22: (1-gamma1*gamma2)*_b[eq7:_cons]) ///
 (a1_cons: _b[eq1:_cons]-gamma1*_b[eq5:_cons]) (a2_cons: _b[eq5:_cons]-gamma1*_b[eq1:_cons]) ///
 (a1gdp: _b[eq1:gdp]-gamma1*_b[eq5:gdp]) (a1gini: _b[eq1:gini]-gamma1*_b[eq5:gini]) ///
 (a1jobs: _b[eq1:jobs]-gamma1*_b[eq5:jobs]) (a1unemp: _b[eq1:unemp]-gamma1*_b[eq5:unemp]) ///
 (a1poprural: _b[eq1:poprural]-gamma1*_b[eq5:poprural]) (a1popmen: _b[eq1:popmen]-gamma1*_b[eq5:popmen]) ///
 (a2gdp: _b[eq5:gdp]-gamma2*_b[eq1:gdp]) (a2gini: _b[eq5:gini]-gamma2*_b[eq1:gini]) ///
 (a2jobs: _b[eq5:jobs]-gamma2*_b[eq1:jobs]) (a2unemp: _b[eq5:unemp]-gamma2*_b[eq1:unemp]) ///
 (a2poprural: _b[eq5:poprural]-gamma2*_b[eq1:poprural]) (a2popmen: _b[eq5:popmen]-gamma2*_b[eq1:popmen]) ///
 (sigma1: exp(_b[eq4:_cons])) (sigma2: exp(_b[eq8:_cons])*(sqrt(6)/_pi)) ///
 (theta: exp(_b[eq9:_cons])/(exp(_b[eq9:_cons])+2)), level(95)

matrix list reduced

matrix b1 = (reduced[1,8..9],reduced[1,1..7])'
matrix b2 = (reduced[1,18..19],reduced[1,11..17])' 
gen cte = 1
mkmat z1 z2 gdp gini jobs unemp poprural popmen cte, matrix(X)
matrix xb1 = X*b1
matrix xb2 = X*b2
svmat double xb1, name(y1star)
svmat double xb2, name(y2star)
rename y1star1 y1star
rename y2star1 y2star
replace y2star = y2star-digamma(1)

/*
tw (scatter y2 y1 if y1star > -.2 & y1 < .8 & y2star > -1.5 & y2 < 2, mcolor(black) msize(small) msymbol(x)), ///
 ytitle(y{subscript:2}) xtitle(y{subscript:1}) ysize(7) xsize(10) ysc(r(-1.5 2)) ylabel(-1.5(.5)2, nogrid) ///
 xsc(r(-.2 .8)) xlabel(-.2(.2).8) graphregion(color(white)) legend(off) xline(0, lpattern(dash) lcolor(gray)) yline(0, lpattern(dash) lcolor(gray))
graph export "C:\Users\Usuario\Meu Drive\artigo_policing\tex\overleaf\s1.png", as(png) replace

tw (scatter y2 y1 if y1> 0 & y1star > -.2 & y1 < .45 & y2 > 0 & y2star > -1.5 & y2 < 1.5, mcolor(black) msize(small) msymbol(x)) ///
 (scatter y2star y1star if y1star > -.2 & y1 < .8 & y2star > -1.5 & y2 < 2, mcolor(black) msize(small) msymbol(x)), ///
 ytitle(predicted y*{subscript:2}) xtitle(predicted y*{subscript:1}) ysize(7) xsize(10) ysc(r(-1.5 2)) ylabel(-1.5(.5)2, nogrid) ///
 xsc(r(-.2 .8)) xlabel(-.2(.2).8) graphregion(color(white)) legend(off) xline(.2, lpattern(dash) lcolor(gray)) yline(0.2, lpattern(dash) lcolor(gray))
graph export "C:\Users\Usuario\Meu Drive\artigo_policing\tex\overleaf\s2.png", as(png) replace
*/

gen mg = gamma1*normal(((y1-y1star)/sigma1-_b[eq9:_cons]*invnormal(exp(-exp(-(y2-y2star)/sigma2))))/sqrt(1-_b[eq9:_cons]^2))
gen elast = mg*y2/y1

tw (kdensity mg) (kdensity elast) if y1>0 & y2 > 0

tw (kdensity mg, sort lcolor(red) lpattern(solid)) ///
 (kdensity elast, sort lcolor(blue) lpattern(solid)) if y1 > 0 & y2 > 0 & elast > -.6, ///
 ytitle(density) xla(-.6(.1)0, axis(1)) ysize(7) xsize(10) graphregion(color(white)) ///
 legend(cols(2) label(1 "marginal effect") label(2 "elasticity")) xtitle("") ylabel(0(2)10, nogrid)
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\elast.eps", as(eps) replace

tw (scatter elast y2, mcolor(blue) msize(small) msymbol(o) ) ///
 (lfit elast y2, lcolor(red) lpattern(solid)) if elast > -.6 & y2 < 1, ///
 ytitle(elasticity) xtitle(y{subscript:2}) ysize(7) xsize(10) ysc(r(-.6 0)) ylabel(-.6(.1)0, nogrid) ///2
 xsc(r(0 1)) xlabel(0(.1)1) graphregion(color(white)) legend(off) 
graph export "E:\OneDrive - UFSC\artigo_policing\tex\overleaf\elast2.eps", as(eps) replace