Introdução à Econometria de Cópulas
O material abaixo traz um mini-curso de “Introdução à Econometria de Cópulas” oferecido para os alunos do CAEN/UFC e da PUC-RS no 2º semestre de 2022, baseado no artigo “Copula econometrics to simulate effects of private policing on crime”, disponível em https://ideas.repec.org/a/ebl/ecbull/eb-20-01192.html, e em outras pesquisas recentes.
Códigos de Stata usados no curso:
cls
clear all
import delimited "C:\Users\Usuario\Meu Drive\curso_copula\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 = 100*1000*hom/pop
gen y2 = 100*psec/police
gen y1_positive = y1 > 0
gen y2_positive = y2 > 0
tab y1_positive y2_positive
su y1 if y1 > 0
su y2 if y2 > 0
tw (scatter y2 y1, mcolor(gray) msize(small) msymbol(o)), ///
ytitle(y2) xtitle(y1) xsize(5) ysize(3) ysc(r(0 375)) ylabel(0(75)375, nogrid) ///
xsc(r(0 125)) xlabel(0(25)125) graphregion(color(white)) legend(off)
su gdp, meanonly
replace gdp = (gdp-r(min))/(r(max)-r(min))
foreach var of varlist gini jobs unemp poprural popmen {
replace `var' = `var'/100
}
replace y1 = y1/100
replace y2 = y2/100
gen z1 = cond(suicides==0,0,fsuicides/suicides)
gen z2 = prison
ktau y1 y2 z1 z2, stats(taub p)
global X gdp gini jobs unemp poprural popmen
sum y1 y2 z1 z2 $X
ivregress 2sls y1 $X (y2 = z2) if y1 > 0 & y2 > 0
scalar alpha1_hat_2sls = e(b)[1,1]
ivregress 2sls y2 $X (y1 = z1) if y1 > 0 & y2 > 0
scalar alpha2_hat_2sls = e(b)[1,1]
reg3 (y1 y2 z1 $X) (y2 y1 z2 $X) if y1 > 0 & y2 > 0
mat list e(b)
scalar alpha1_hat_reg3 = e(b)[1,1]
scalar alpha2_hat_reg3 = e(b)[1,10]
sureg (y1 z1 z2 $X) (y2 z1 z2 $X) if y1 > 0 & y2 > 0
mat list e(b)
scalar alpha1_hat_sur = e(b)[1,2]/e(b)[1,11]
scalar alpha2_hat_sur = e(b)[1,10]/e(b)[1,1]
disp alpha1_hat_2sls
disp alpha1_hat_sur
disp alpha1_hat_reg3
disp alpha2_hat_2sls
disp alpha2_hat_sur
disp alpha2_hat_reg3
program define fiml_gaussian
args lnf xb1 xb2 s1 s2 r
tempvar v1 v2 rho
quietly {
local y1 "$ML_y1"
local y2 "$ML_y2"
gen double `v1' = (`y1'-`xb1')/exp(`s1')
gen double `v2' = (`y2'-`xb2')/exp(`s2')
gen double `rho' = (2/_pi)*atan(`r')
replace `lnf' = -`s1'-`s2'-ln(2*_pi)-.5*ln(1-`rho'^2)-(.5/(1-`rho'^2))*(`v1'^2-2*`rho'*`v1'*`v2'+`v2'^2)
}
end
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml check
constraint 1 _b[eq5:_cons] = 0
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0, constraint(1)
ml max, difficult iterate(50)
matrix start = e(b)
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml init start
ml max, difficult iterate(50)
cls
clear all
set obs 10000
matrix g = J(sqrt(_N),1,0)
local b = 4
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
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_gI = exp(-exp(-vv1))
gen g1_gI = exp(-vv1-exp(-vv1))/b
gen vv2 = (v2-a)/b
gen G2_gI = exp(-exp(-vv2))
gen g2_gI = exp(-vv2-exp(-vv2))/b
replace vv1 = (v1+a)/b
replace vv2 = (v2+a)/b
gen G1_gII = 1-exp(-exp(vv1))
gen g1_gII = exp(vv1-exp(vv1))/b
gen G2_gII = 1-exp(-exp(vv2))
gen g2_gII = exp(vv2-exp(vv2))/b
tw (line g1_gI v1, sort lcolor(gray) lpattern(shortdash)) ///
(line g1_gauss v1, sort lcolor(black) lpattern(solid)) ///
(line g1_gII v1, sort lcolor(black) lpattern(dash)), ///
ytitle(density) xtitle(variable) xla(-4(1)4, axis(1)) xsize(6) ysize(4) scheme(s1mono) ///
legend(cols(3) label(1 "Gumbel I") label(2 "Normal") label(3 "Gumbel II"))
scalar theta = .707 /* tau = .5 */
gen C12 = (exp(-(invnormal(G1_gI)^2-2*theta*invnormal(G1_gI)*invnormal(G2_gauss)+invnormal(G2_gauss)^2)/(2*(1-theta^2)))/(2*_pi*sqrt(1-theta^2)))/ ///
(normalden(invnormal(G1_gI))*normalden(invnormal(G2_gauss)))
gen f12 = g1_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
replace C12 = (exp(-(invnormal(G1_gI)^2-2*theta*invnormal(G1_gI)*invnormal(G2_gI)+invnormal(G2_gI)^2)/(2*(1-theta^2)))/(2*_pi*sqrt(1-theta^2)))/ ///
(normalden(invnormal(G1_gI))*normalden(invnormal(G2_gI)))
replace f12 = g1_gI*g2_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
replace C12 = (exp(-(invnormal(G1_gII)^2-2*theta*invnormal(G1_gII)*invnormal(G2_gI)+invnormal(G2_gI)^2)/(2*(1-theta^2)))/(2*_pi*sqrt(1-theta^2)))/ ///
(normalden(invnormal(G1_gII))*normalden(invnormal(G2_gI)))
replace f12 = g1_gII*g2_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
scalar theta = 2 /* tau = .5 */
gen C = (G1_gI^(-theta)+G2_gauss^(-theta)-1)^(-1/theta)
replace C12 = ((1+theta)*(C^(1+2*theta)))/((G1_gI*G2_gauss)^(1+theta))
replace f12 = g1_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
replace C = (G1_gI^(-theta)+G2_gI^(-theta)-1)^(-1/theta)
replace C12 = ((1+theta)*(C^(1+2*theta)))/((G1_gI*G2_gI)^(1+theta))
replace f12 = g1_gI*g2_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
replace C = (G1_gII^(-theta)+G2_gI^(-theta)-1)^(-1/theta)
replace C12 = ((1+theta)*(C^(1+2*theta)))/((G1_gII*G2_gI)^(1+theta))
replace f12 = g1_gII*g2_gI*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(v1) ytitle(v2) ///
xla(-4(1)4, axis(1)) yla(-4(1)4, axis(1)) ysize(7) xsize(10) scheme(s1mono)
cls
clear all
import delimited "C:\Users\Usuario\Meu Drive\curso_copula\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 %"
su gdp, meanonly
replace gdp = (gdp-r(min))/(r(max)-r(min))
foreach var of varlist gini jobs unemp poprural popmen {
replace `var' = `var'/100
}
gen y1 = 1000*hom/pop
gen y2 = psec/police
gen z1 = cond(suicides==0,0,fsuicides/suicides)
gen z2 = prison
global X gdp gini jobs unemp poprural popmen
program define fiml_gaussian
args lnf xb1 xb2 s1 s2 r
tempvar v1 v2 rho
quietly {
local y1 "$ML_y1"
local y2 "$ML_y2"
gen double `v1' = (`y1'-`xb1')/exp(`s1')
gen double `v2' = (`y2'-`xb2')/exp(`s2')
gen double `rho' = (2/_pi)*atan(`r')
replace `lnf' = -`s1'-`s2'-ln(2*_pi)-.5*ln(1-`rho'^2)-(.5/(1-`rho'^2))*(`v1'^2-2*`rho'*`v1'*`v2'+`v2'^2)
}
end
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml check
constraint 1 _b[eq5:_cons] = 0
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0, constraint(1)
ml max, difficult iterate(50)
matrix start = e(b)
ml model lf fiml_gaussian (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml init start
ml max, difficult iterate(50)
estimates store gauss_n_n
program define fiml_clayton
args lnf xb1 xb2 s1 s2 r
tempvar v1 v2 theta G1 g1 G2 g2 D C12
quietly {
local y1 "$ML_y1"
local y2 "$ML_y2"
gen double `v1' = (`y1'-`xb1')/exp(`s1')
gen double `v2' = (`y2'-`xb2')/exp(`s2')
gen double `theta' = exp(`r')
gen double `G1' = exp(-exp(-`v1'))
gen double `g1' = exp(-`v1'-exp(-`v1'))
gen double `G2' = exp(-exp(-`v2'))
gen double `g2' = exp(-`v2'-exp(-`v2'))
gen double `D' = (`G1'^(-`theta')+`G2'^(-`theta')-1)^(-1/`theta')
gen double `C12' = ((1+`theta')/((`G1'*`G2')^(1+`theta')))*(`D'^(1+2*`theta'))
replace `lnf' = ln(max(1e-20,(`g1'/exp(`s1'))*(`g2'/exp(`s2'))*`C12'))
}
end
ml model lf fiml_clayton (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml check
constraint 1 _b[eq5:_cons] = 0
ml model lf fiml_clayton (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0, constraint(1)
ml max, difficult iterate(50)
matrix start = e(b)
ml model lf fiml_clayton (y1 y2 = z1 z2 $X) (z1 z2 $X) () () () if y1 > 0 & y2 > 0
ml init start
ml max, difficult iterate(50)
estimates store clayton_gI_gI
esttab gauss_n_n clayton_gI_gI, se stats(ll)
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 = 1
gen k = 1
/*
*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 y2star y1star if y1star > -.2 & y1 < .8 & y2star > -1.5 & y2 < 2, mcolor(black) msize(tiny) msymbol(x)) ///
(scatter y2 y1 if y1star > -.2 & y1 < .8 & y2star > -1.5 & y2 < 2, mcolor(gray) msize(small) msymbol(dh)), ///
ytitle(y2*) xtitle(y1*) xsize(10) ysize(6) ysc(r(-1.5 2)) ylabel(-1.5(.5)2, nogrid) ///
xsc(r(-.2 .8)) xlabel(-.2(.2).8) graphregion(color(white)) legend(off)
graph export "C:\Users\Usuario\Meu Drive\artigo_policing\tex\scatter2.png", as(png) replace
***********************************
tobit y1 $X z1 z2, ll(0)
matrix b = e(b)
matrix start1 = (b[1..1,1..6],b[1,9],b[1,7],b[1,8])
scalar sd1 = ln(sqrt(b[1,10]))
tobit y2 $X z1 z2, ll(0)
matrix b = e(b)
matrix start2 = (b[1..1,1..6],b[1,9],b[1,7],b[1,8])
scalar sd2 = ln(sqrt(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