Julia 101: uma Abordagem Regional
Primeiramente, vamos instalar e ler as bibliotecas necessárias para o workshop
using Pkg # pacote para instalar novos pacotes
using LinearAlgebra # pacote de Álgebra Linear
using SpatialDependence # pacote de análise de dados espaciais
using Optim # pacote com funções básicas de otimização
using JuMP,Ipopt # pacote com algoritmos de otimização com restrições
using Random # para criar vetores, matrizes e arranjos aleatórios
using Shapefile # para trabalhar com shapefiles
using Plots # para plot de shapefiles simples
using Jenks # para quebra natural de classes em jenks
using StableRNGs # para travar o gerador de números aleatórios na análise espacial
using DataFrames # para usar base de dados em formatos retangulares
using ForwardDiff # para calcular segundas derivadas
using ProgressBars # para barras de progresso
using PrettyTables # para criar tabelas apresentáveis
using Statistics # para funções estatísticas
Como instalamos um novo pacote na Julia? Usando o comando Pkg.add().
Primeiramente, criemos uma matriz de números aleatórios, distribuídos conforme U[0,1].
Random.seed!(12345)
A=rand(2,2)
2×2 Matrix{Float64}:
0.944791 0.339612
0.866895 0.136117
Alguém se lembra da fórmula para criar números distribuídos uniformemente a partir de uma U[0,1]?
B=rand(1.0:10.0,2,2)
C=rand(1:10,2,2)
B,C
([4.0 7.0; 3.0 8.0], [4 1; 6 9])
Por que B é do tipo float, enquanto C é do tipo int? Aqui a typagem forte da Julia fica mais explícita. Observe que a matriz A anterior também é float, pois ela puxa a distribuição uniforme [0.0,1.0]. Vamos agora realizar duas operações que são muito comuns na Econometria: inversão de matrizes e encontrar seus auto-valores.
inv(C)
2×2 Matrix{Float64}:
0.3 -0.0333333
-0.2 0.133333
eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
-0.1362265186681017
1.217133601568467
vectors:
2×2 Matrix{Float64}:
-0.299717 0.780136
0.954028 0.62561
Vamos agora explorar um algoritmo otimização não condicional na Julia. Considere a segunte função: $y=-(x_{1}^2+x_{2}^2+x_{3}^{2})$
function f(x)
return +(x[1]^2 + x[2]^2 + x[3]^2)
end
# Usando agora alguns pontos iniciais|
x0 = [1.0, 1.0, 1.0]
# Realizando agora a otimização
resultado = optimize(f, x0)
resultado
* Status: success
* Candidate solution
Final objective value: 2.833427e-09
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 67
f(x) calls: 123
Para ver os pontos que minimizam a função, basta usar o qualificador minimizer
resultado.minimizer
3-element Vector{Float64}:
-5.0823878310249317e-5
9.262030970117879e-6
-1.2828696727829803e-5
Para ver o valor máximo que sua função adquire, baste usar o qualificador minimum
resultado.minimum
2.8334272839210848e-9
Um detalhe importante é que os algoritmos de otimização do pacote Optim minimizam sua função a ser otimizada. Então, é importante garantir que ela tem solução e que essa função entra no processo de otimização de forma correta. Como uma otimização não condicionada funcionaria. Usemos o caso de uma função de utilidade do tipo Cobb-Douglas simples para entender a dinâmica do código. Aqui, vamos usar o pacote JUMP que permite otimizações com restrições
# Criando o modelo (vazio)
modelo = Model(Ipopt.Optimizer)
# Definir as variáveis - não-negativas
@variable(modelo, x1 >= 0)
@variable(modelo, x2 >= 0)
# Define the objective function to maximize
@objective(modelo, Max, (x1^2 + x2^2))
# Add linear constraints
@constraint(modelo, 2*x1 + 5*x2 <= 100.0)
# Optimize the model
optimize!(modelo)
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 2
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 2
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.9999960e-04 0.00e+00 7.29e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 2.5067536e-02 0.00e+00 2.18e-01 -1.0 7.13e-01 - 8.20e-01 1.00e+00f 1
2 2.6040386e-02 0.00e+00 6.02e-01 -1.0 6.06e-03 2.0 9.95e-01 1.00e+00f 1
3 6.7549927e-02 0.00e+00 3.92e-01 -1.7 4.75e-01 - 6.92e-01 1.00e+00f 1
4 6.9927597e-02 0.00e+00 2.84e-01 -1.7 8.51e-03 1.5 1.00e+00 1.00e+00f 1
5 7.9112567e-02 0.00e+00 3.21e-01 -1.7 2.89e-02 1.0 1.00e+00 1.00e+00f 1
6 1.5617712e-01 0.00e+00 5.58e-01 -1.7 1.51e-01 0.6 1.00e+00 1.00e+00f 1
7 2.0927484e-01 0.00e+00 6.91e-01 -1.7 6.99e-02 1.0 1.00e+00 1.00e+00f 1
8 4.5448940e-01 0.00e+00 1.02e+00 -1.7 3.95e-01 0.5 1.00e+00 5.68e-01f 1
9 9.7151939e-01 0.00e+00 8.28e-01 -1.7 7.55e-01 0.0 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.4369612e+00 0.00e+00 8.93e-01 -1.7 3.05e-01 0.5 1.00e+00 1.00e+00f 1
11 4.1529868e+00 0.00e+00 1.63e+00 -1.7 1.67e+00 -0.0 1.00e+00 1.00e+00f 1
12 5.8013089e+00 0.00e+00 1.93e+00 -1.7 7.42e-01 0.4 1.00e+00 1.00e+00f 1
13 2.0054598e+01 0.00e+00 3.58e+00 -1.7 4.13e+00 -0.1 1.00e+00 1.00e+00f 1
14 2.9306227e+01 0.00e+00 4.33e+00 -1.7 1.87e+00 0.4 1.00e+00 1.00e+00f 1
15 1.0431182e+02 0.00e+00 8.17e+00 -1.7 1.17e+01 -0.1 8.15e-01 8.21e-01f 1
16 1.6079098e+02 0.00e+00 1.01e+01 -1.7 4.94e+00 0.3 1.00e+00 1.00e+00f 1
17 5.6248055e+02 0.00e+00 1.90e+01 -1.7 3.56e+01 -0.2 6.21e-01 6.20e-01f 1
18 9.2200033e+02 0.00e+00 2.43e+01 -1.7 1.33e+01 0.3 1.00e+00 1.00e+00f 1
19 2.4803885e+03 0.00e+00 3.98e+01 -1.7 1.16e+02 -0.2 4.65e-01 3.34e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 2.4997886e+03 0.00e+00 6.43e+01 -1.7 3.25e+01 0.2 1.00e+00 1.20e-02f 1
21 2.4999829e+03 0.00e+00 5.29e+01 -1.7 1.24e+00 - 1.00e+00 3.14e-03f 1
22 2.4999603e+03 0.00e+00 2.00e-07 -1.7 3.56e-04 - 1.00e+00 1.00e+00f 1
23 2.4999998e+03 0.00e+00 1.50e-09 -3.8 3.94e-04 - 1.00e+00 1.00e+00f 1
24 2.5000000e+03 0.00e+00 1.84e-11 -5.7 2.97e-06 - 1.00e+00 1.00e+00f 1
25 2.5000001e+03 0.00e+00 5.68e-14 -8.6 3.68e-08 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 25
(scaled) (unscaled)
Objective...............: -2.5000000524949887e+03 2.5000000524949887e+03
Dual infeasibility......: 5.6843418860808015e-14 5.6843418860808015e-14
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 9.9899763890181942e-09 9.9899763890181942e-09
Complementarity.........: 2.5060842553385023e-09 2.5060842553385023e-09
Overall NLP error.......: 2.5060842288585504e-09 2.5060842553385023e-09
Number of objective function evaluations = 26
Number of objective gradient evaluations = 26
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 26
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 0.138
EXIT: Optimal Solution Found.
Vamos agora extrair nossos objetos de interesse. O resultado indica que o consumidor só consome o bem 1, mais barato e que lhe produz tanta utilidade quanto o bem 1
value(x1),value(x2),objective_value(modelo)
(50.00000052494988, -9.989976389018194e-9, 2500.0000524949887)
Em algumas situações, é possível transferir um problema de otimização com restrições em algum outro problema de otimização. Logo, não é sempre necessário usar o pacote JUMP, bastando o próprio Optim. Ao mesmo tempo, há situações de otimizações numéricas mais complexas que o pacote JUMP não é suficiente, contudo há uma grande quantidade de pacotes para essa finalidade no escossistema da Julia
Agora, vamos usar o pacote SpatialDependece para analisar correlação espacial para algumas variáveis relacionadas a Belo Horizonte-MG. Criaremos matrizes espaciais e plotaremos alguns mapas neste exercício. Inicialmente, leiamos o shapefile necessário.
shape = Shapefile.Table("/Users/alanleal/Library/CloudStorage/OneDrive-Personal/Aulas/Julia 101 uma abordagem regional/Shapefiles/POP_DOMIC_BAIRRO_2010/POP_DOMIC_BAIRRO_2010.shp")
plot(shape)
Vamos agora criar uma matriz de pesos espacial
W = polyneigh(shape.geometry)
W
Spatial Weights Observations: 487 Transformation: row Minimum nunmber of neighbors: 0 Maximum nunmber of neighbors: 17 Average number of neighbors: 4.7803 Median number of neighbors: 5.0 Islands (isloated): 57 Density: 0.9816%
Calculando agora o I de Moran
moran(shape.POPULACAO, W, permutations = 9999, rng = StableRNG(1234567))
Moran's I test of Global Spatial Autocorrelation -------------------------------------------- Moran's I: 0.0603264 Expectation: -0.0020576 Randomization test with 9999 permutations. Mean: -0.0025178 Std Error: 0.031551 zscore: 1.9918235 p-value: 0.0282
Plotando agora o I de Moran
plot(shape.POPULACAO, W, xlabel = "População")
Agora, calculando e plotando o LISA
lmpop = localmoran(shape.POPULACAO, W, permutations = 9999, rng = StableRNG(1234567))
plot(shape, lmpop, sig = 0.05, adjust = :fdr)
Mais detalhes, inclusive outros índices de (auto-)correlação espacial e diferentes matrizes espaciais, podem ser conferidos no seguinte link: https://javierbarbero.github.io/SpatialDependence.jl/stable/.
Nosso objetivo aqui é implementar um simples estimador da Econometria Espacial, qual seja o SAR. Vamos implementar três coisas inicialmente; (i) estimação pontual dos coeficientes de interesse; (ii) estimação intervalar dos coeficientes de interesse; (iii) efeitos diretos, indiretos e totais das observações
4.1. Leitura do Shapefile - Casos de Zika no Ceará e Mapas
Este shapefile foi obtido no seguinte link:https://geodacenter.github.io/data-and-lab/Ceara-Zika/
# Origem da base: https://geodacenter.github.io/data-and-lab/
ceara_zika=Shapefile.Table("/Users/alanleal/Library/CloudStorage/OneDrive-Personal/Aulas/Julia 101 uma abordagem regional/Shapefiles/ceara/ceara.shp")
plot(ceara_zika)
Num segundo momento, vamos plotar um mapa de cloropeth mais tradicional, qual seja um mapa de quantiles. Para tanto, criamos uma nova coluna identificando a qual quantil cada município do Ceará pertence.
df=DataFrame(ceara_zika)
quantiles = quantile(df.ln_gdp, 0:0.2:1)
function classify_quantiles(data, quantiles)
return [findfirst(q -> x <= q, quantiles) for x in data]
end
df[!,:quantile_class] = classify_quantiles(df[!,:ln_gdp], quantiles)
plot(df.geometry, fill_z=reshape(df.quantile_class,1,nrow(df)), cbar=true, axis=false, ticks=false, size=(800, 600))
Agora, plotando o mesmo mapa cloropeth com quebras naturais (jenks)
# Transformando um vetor de missings e float64 em float64 apenas
mean_value = mean(skipmissing(df.ln_gdp))
df.ln_gdp = replace(df.ln_gdp, missing => mean_value)
plt_bounds=JenksClassification(5,df.ln_gdp).bounds
function classify_jenks(data, jenks)
return [findfirst(q -> x <= q, jenks) for x in data]
end
df[!,:jenks_class] = classify_jenks(df[!,:ln_gdp], plt_bounds)
plot(df.geometry, fill_z=reshape(df.jenks_class,1,nrow(df)), cbar=true, axis=false, ticks=false, size=(800, 600))
100%: ARE=0.10351076, finished in 0.01s.
4.2. Criação de objetos de interesse
Inicialmente, criemos uma matriz normalizada de continguidade
W=polyneigh(ceara_zika.geometry) # matriz de contiguidade do fomato queen-default
# W = polyneigh(ceara_zika.geometry, criterion = :Rook) # matriz de contiguidade do fomato Rook
# W = dnearneigh(ceara_zika.geometry, threshold = 4.0) # matriz de inverso da distância para shapefile de PONTOS
# W = knearneigh(boston.geometry, k = 5) # matriz de k=5 vizinhos mais próximos para shapefile de PONTOS
W=Matrix(W) # transforma a matriz espacial em uma matrix usual
W
184×184 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 … 0.25 0.25 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.142857 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Selecionando as variáveis para análise
n=184 # número de observações
y = ceara_zika.inc_zik_3q # jeito de selecionar uma variável do shapefile ceara_zika
df = DataFrame(
constant=ones(n),
ln_gdp = ceara_zika.ln_gdp,
ln_pop=ceara_zika.ln_pop,
mobility = ceara_zika.mobility,
environ = ceara_zika.environ,
sanitation = ceara_zika.sanitation
)
X = Matrix(select(df, [:constant,:ln_gdp,:ln_pop,:mobility, :environ, :sanitation]))
X
184×6 Matrix{Union{Missing, Float64}}:
1.0 4.55599 4.02102 0.958 0.835 0.538
1.0 4.83451 4.18577 0.941 0.644 0.728
1.0 5.49065 4.76005 0.971 0.973 0.535
1.0 5.31389 4.70893 0.975 0.885 0.687
1.0 4.73446 4.2096 0.99 0.96 0.68
1.0 4.57723 4.03226 1.0 0.932 0.798
1.0 4.3909 3.83607 0.968 0.538 0.644
1.0 4.88703 4.21376 0.972 0.905 0.573
1.0 5.26288 4.59364 0.933 0.877 0.524
1.0 4.45245 3.8441 0.956 0.926 0.752
⋮ ⋮
1.0 4.71922 4.1586 0.95 0.859 0.564
1.0 5.29099 4.50225 0.981 0.788 0.607
1.0 4.46604 3.87766 0.976 0.837 0.625
1.0 4.85375 4.2742 0.933 0.8 0.589
1.0 5.14761 4.2959 0.943 0.906 0.612
1.0 4.68884 4.11002 0.978 0.938 0.659
1.0 4.99156 4.24534 0.98 0.917 0.586
1.0 5.20945 4.58472 0.956 0.823 0.679
1.0 5.3319 4.74001 0.973 0.77 0.557
Agora, criemos a função de log-verossimilhança do modelo SAR. Ele é dado por:
$l(\theta|y,X,W) \propto -(n/2)*log(\varepsilon'\varepsilon)+log(|I-\rho W|)$, com $\varepsilon=(I-\rho W)(y-X\beta)$. Com base nisso, criamos a seguinte função que avalia a log-verossimilhança do modelo SAR
function sar_likelihood(params)
ρ, β = params[1], params[2:end]
ε = (I(n) - ρ*W)*y-X*β
loglik = ((-n/2)*log(ε'*ε))+log(det(I(n)-ρ*W))
return -loglik # Negativo da log-verossimilhança para de fato maximizar essa função
end
sar_likelihood (generic function with 1 method)
Vejamos agora se a função de verossimilhança avalia os dados corretamente
sar_likelihood([0,0,0,0,0,0,0])
960.775240104947
Agora, vamos considerar os dados de Zika no Ceará para estimar um simples modelo espacial das variáveis que impactam os casos municipais de Zika no terceiro trimestre de 2019. Lembre-se que enquanto $\beta$ é livre, o parâmetro $\rho$ é limitado entre [-1,1], ou ainda entre $[min[-1/(min_{i}\{\lambda_{i}\})],1]$, com $min_{i}\{\lambda_{i}\}$ sendo o menor autovalor da matriz W de pesos espaciais. Ademais, no comando a seguir, temos quatro componentes, respectivamente: (i) parâmetros iniciais, para começo da busca; (ii) lower bounds, que indicam os valores mínimos assumidos por cada parâmetro; (iii) upper bounds, que indicam os valores máximos para cada parâmetro; (iv) result, que invoca a função optimize, que minimiza, do pacote Optim
initial_params = [0.5,0,0,0,0,0,0] # Initial values for ρ e β
lower_bounds = [-1, -Inf,-Inf,-Inf,-Inf,-Inf,-Inf]
upper_bounds = [1, Inf,Inf,Inf,Inf,Inf,Inf]
result = optimize(sar_likelihood,lower_bounds, upper_bounds,initial_params,Fminbox())
* Status: success
* Candidate solution
Final objective value: 9.565912e+02
* Found with
Algorithm: Fminbox with L-BFGS
* Convergence measures
|x - x'| = 8.67e-08 ≰ 0.0e+00
|x - x'|/|x'| = 3.75e-09 ≰ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 9.39e-09 ≤ 1.0e-08
* Work counters
Seconds run: 1 (vs limit Inf)
Iterations: 3
f(x) calls: 133
∇f(x) calls: 133
Acessemos agora nossos coeficientes β, com a opção minimizer
β = result.minimizer
β
7-element Vector{Float64}:
0.23790174678155787
16.62706670088536
-2.59958822275927
1.2283227903363063
-14.589927485596965
5.229097057229507
3.2474406362839607
Num segundo momento, calculemos agora a matriz de variância e covarância do SAR. Ela é dada por: $\hat{\omega}=\sigma^{2}(I-\rho W)^{-1}(I-\rho W')^{-1}$. Calculemos $\sigma^{2}$ inicialmente
haty=β[1]*W*y+X*β[2:end] #y estimado
sigma2=sum((y-haty).^(2))/(184-6) # erro amostral dos resíduos
sigma2
182.1539950249789
Por fim, calculando a matriz de variâncias e covariâncias dos parâmetros. Esse cálculo usa o cálculo de uma segunda derivada numericamente. Para detalhamento teórico, veja: https://www.researchgate.net/publication/229020633_Lecture_1_Maximum_likelihood_estimation_of_spatial_regression_models
σk = inv(ForwardDiff.hessian(sar_likelihood, β))
sqrt.(diag(σk))
7-element Vector{Float64}:
0.11315758632823769
39.88016539980641
6.923082060731218
8.970209233447832
35.78905507302816
11.160015889983043
10.118056357643486
Vamos agora agora colocar os coeficientes com seus respectivos desvios-padrão numa única matriz
hcat(β,sqrt.(diag(σk)))
7×2 Matrix{Float64}:
0.237902 0.113158
16.6271 39.8802
-2.59959 6.92308
1.22832 8.97021
-14.5899 35.7891
5.2291 11.16
3.24744 10.1181
4.3. Efeitos Diretos e Indiretos: Estimação Pontual
ψ=inv(I(n)-β[1]*W)
eft_direto=zeros(length(β)-1)
eft_indireto=zeros(length(β)-1)
eft_total=zeros(length(β)-1)
for i=2:length(β)
eft_direto[i-1]=(1/n)*tr(inv(I(n)-β[1]*W))β[i]
eft_total[i-1]=(1/n)*sum(sum(ψ[k,j]*β[i] for j=1:n) for k=1:n)
eft_indireto[i-1]=eft_total[i-1]-eft_direto[i-1]
end
Vejamos agora esse três vetores:
Efeitos=hcat(eft_direto,eft_indireto,eft_total)
Efeitos
6×3 Matrix{Float64}:
16.8187 4.99875 21.8175
-2.62956 -0.781538 -3.41109
1.24248 0.369282 1.61176
-14.7581 -4.3863 -19.1444
5.28938 1.57207 6.86145
3.28488 0.976308 4.26118
4.4. Efeitos Diretos e Indiretos: Estimação Intervalar via Bootstrap
Agora, defininamos algumas funções de interesse antes. A função a seguir calcula os efeitos diretos, indiretos e totais, considerando um vetor de parâmetros $\beta_{1}$. Ela será usada continuamente na nossa função de bootstrap.
function efeitos(β1)
ψ=inv(I(n)-β1[1]*W)
eft_direto=zeros(length(β1)-1)
eft_indireto=zeros(length(β1)-1)
eft_total=zeros(length(β1)-1)
for i=2:length(β1)
eft_direto[i-1]=(1/n)*tr(inv(I(n)-β1[1]*W))*β1[i]
eft_total[i-1]=(1/n)*sum(sum(ψ[k,j]*β1[i] for j=1:n) for k=1:n)
eft_indireto[i-1]=eft_total[i-1]-eft_direto[i-1]
end
return hcat(eft_direto,eft_indireto,eft_total)
end
efeitos (generic function with 1 method)
Agora, definamos a função de Bootstrap
# Inicialmente, defina o número de bootstraps.
n_bootstrap = 50
# Criando objetos vazios para acumular com nossas estimativas pontuais
bootstrap_direct = []
bootstrap_indirect = []
bootstrap_total = []
y2=y
X2=X
for i in ProgressBar(1:n_bootstrap)
# Reamostragem com reposição
resample_indices = rand(1:length(y), length(y))
Y_resample = y2[resample_indices]
X_resample = X2[resample_indices, :]
# Otimiza a função SAR em cada iteração do bootstrap
X=X_resample
y=Y_resample
β1 = optimize(sar_likelihood,lower_bounds, upper_bounds,initial_params,Fminbox()).minimizer
# Calcula os efeitos com base na função que definimos anteriormente
efeitos_df=efeitos(β1)
push!(bootstrap_direct, efeitos_df[:,1])
push!(bootstrap_indirect, efeitos_df[:,2])
push!(bootstrap_total, efeitos_df[:,3])
end
# Por fim, captura o desvio-padrão calculado via bootstrap
std_direct = std(bootstrap_direct)
std_indirect = std(bootstrap_indirect)
std_total = std(bootstrap_total)
0.0%┣ ┫ 0/50 [00:01<00:-38, -1s/it] 2.0%┣▉ ┫ 1/50 [00:02<Inf:Inf, InfGs/it] 4.0%┣██ ┫ 2/50 [00:03<02:33, 3s/it] 6.0%┣██▉ ┫ 3/50 [00:04<01:44, 2s/it] 8.0%┣███▉ ┫ 4/50 [00:06<01:30, 2s/it] 10.0%┣████▊ ┫ 5/50 [00:07<01:19, 2s/it] 12.0%┣█████▋ ┫ 6/50 [00:09<01:19, 2s/it] 14.0%┣██████▋ ┫ 7/50 [00:10<01:11, 2s/it] 16.0%┣███████▌ ┫ 8/50 [00:11<01:06, 2s/it] 18.0%┣████████▌ ┫ 9/50 [00:12<01:04, 2s/it] 20.0%┣█████████▏ ┫ 10/50 [00:13<00:59, 1s/it] 22.0%┣██████████▏ ┫ 11/50 [00:14<00:55, 1s/it] 24.0%┣███████████ ┫ 12/50 [00:15<00:52, 1s/it] 26.0%┣████████████ ┫ 13/50 [00:16<00:49, 1s/it] 28.0%┣████████████▉ ┫ 14/50 [00:17<00:46, 1s/it] 30.0%┣█████████████▉ ┫ 15/50 [00:18<00:45, 1s/it] 32.0%┣██████████████▊ ┫ 16/50 [00:19<00:43, 1s/it] 34.0%┣███████████████▋ ┫ 17/50 [00:20<00:42, 1s/it] 36.0%┣████████████████▋ ┫ 18/50 [00:21<00:40, 1s/it] 38.0%┣█████████████████▌ ┫ 19/50 [00:22<00:38, 1s/it] 40.0%┣██████████████████▍ ┫ 20/50 [00:23<00:37, 1s/it] 42.0%┣███████████████████▎ ┫ 21/50 [00:24<00:36, 1s/it] 44.0%┣████████████████████▎ ┫ 22/50 [00:25<00:34, 1s/it] 46.0%┣█████████████████████▏ ┫ 23/50 [00:26<00:32, 1s/it] 48.0%┣██████████████████████ ┫ 24/50 [00:27<00:31, 1s/it] 50.0%┣███████████████████████ ┫ 25/50 [00:28<00:30, 1s/it] 52.0%┣████████████████████████ ┫ 26/50 [00:29<00:28, 1s/it] 54.0%┣████████████████████████▉ ┫ 27/50 [00:30<00:27, 1s/it] 56.0%┣█████████████████████████▊ ┫ 28/50 [00:32<00:26, 1s/it] 58.0%┣██████████████████████████▊ ┫ 29/50 [00:33<00:25, 1s/it] 60.0%┣███████████████████████████▋ ┫ 30/50 [00:34<00:23, 1s/it] 62.0%┣████████████████████████████▌ ┫ 31/50 [00:34<00:22, 1s/it] 64.0%┣█████████████████████████████▍ ┫ 32/50 [00:35<00:20, 1s/it] 66.0%┣██████████████████████████████▍ ┫ 33/50 [00:36<00:19, 1s/it] 68.0%┣███████████████████████████████▎ ┫ 34/50 [00:37<00:18, 1s/it] 70.0%┣████████████████████████████████▏ ┫ 35/50 [00:38<00:17, 1s/it] 72.0%┣█████████████████████████████████▏ ┫ 36/50 [00:40<00:16, 1s/it] 74.0%┣██████████████████████████████████ ┫ 37/50 [00:41<00:15, 1s/it] 76.0%┣███████████████████████████████████ ┫ 38/50 [00:42<00:14, 1s/it] 78.0%┣███████████████████████████████████▉ ┫ 39/50 [00:43<00:12, 1s/it] 80.0%┣████████████████████████████████████▉ ┫ 40/50 [00:44<00:11, 1s/it] 82.0%┣█████████████████████████████████████▊ ┫ 41/50 [00:45<00:10, 1s/it] 84.0%┣██████████████████████████████████████▋ ┫ 42/50 [00:46<00:09, 1s/it] 86.0%┣███████████████████████████████████████▋ ┫ 43/50 [00:47<00:08, 1s/it] 88.0%┣████████████████████████████████████████▌ ┫ 44/50 [00:48<00:07, 1s/it] 90.0%┣█████████████████████████████████████████▍ ┫ 45/50 [00:49<00:06, 1s/it] 92.0%┣██████████████████████████████████████████▎ ┫ 46/50 [00:50<00:04, 1s/it] 94.0%┣███████████████████████████████████████████▎ ┫ 47/50 [00:51<00:03, 1s/it] 96.0%┣████████████████████████████████████████████▏ ┫ 48/50 [00:53<00:02, 1s/it] 98.0%┣█████████████████████████████████████████████ ┫ 49/50 [00:54<00:01, 1s/it] 100.0%┣█████████████████████████████████████████████┫ 50/50 [00:55<00:00, 1s/it] 100.0%┣█████████████████████████████████████████████┫ 50/50 [00:55<00:00, 1s/it]
6-element Vector{Float64}:
26.662380198845682
2.7265856146621115
3.4543689988386976
28.3518226484923
6.681105575085341
3.5183713295663144
Efeitos Diretos
hcat(Efeitos[:,1],std_direct)
6×2 Matrix{Float64}:
16.8187 26.1925
-2.62956 2.76359
1.24248 3.51387
-14.7581 27.9545
5.28938 6.63234
3.28488 3.50714
Efeitos Indiretos
hcat(Efeitos[:,2],std_indirect)
6×2 Matrix{Float64}:
4.99875 2.47909
-0.781538 0.208176
0.369282 0.163047
-4.3863 2.34135
1.57207 0.578201
0.976308 0.293798
Efeitos Totais
hcat(Efeitos[:,3],std_total)
6×2 Matrix{Float64}:
21.8175 26.6624
-3.41109 2.72659
1.61176 3.45437
-19.1444 28.3518
6.86145 6.68111
4.26118 3.51837
Construamos agora uma tabela com todos esses efeitos
col_names = ["Variável","Efeitos Diretos", "Desvio-Padrão ED", "Efeitos Indiretos", "Desvio-Padrão EI", "Efeitos Totais", "Desvio-Padrão ET"]
df = DataFrame(hcat(["Constant","ln_gdp","ln_pop","mobility","environ","sanitation"],Efeitos[:,1],std_direct,Efeitos[:,2],std_indirect,Efeitos[:,3],std_total), col_names)
pretty_table(df, header = col_names, alignment = :c, formatters = ft_printf("%.2f", 1:6))
┌────────────┬─────────────────┬──────────────────┬───────────────────┬──────────────────┬────────────────┬──────────────────┐ │ Variável │ Efeitos Diretos │ Desvio-Padrão ED │ Efeitos Indiretos │ Desvio-Padrão EI │ Efeitos Totais │ Desvio-Padrão ET │ ├────────────┼─────────────────┼──────────────────┼───────────────────┼──────────────────┼────────────────┼──────────────────┤ │ Constant │ 16.82 │ 26.19 │ 5.00 │ 2.48 │ 21.82 │ 26.6624 │ │ ln_gdp │ -2.63 │ 2.76 │ -0.78 │ 0.21 │ -3.41 │ 2.72659 │ │ ln_pop │ 1.24 │ 3.51 │ 0.37 │ 0.16 │ 1.61 │ 3.45437 │ │ mobility │ -14.76 │ 27.95 │ -4.39 │ 2.34 │ -19.14 │ 28.3518 │ │ environ │ 5.29 │ 6.63 │ 1.57 │ 0.58 │ 6.86 │ 6.68111 │ │ sanitation │ 3.28 │ 3.51 │ 0.98 │ 0.29 │ 4.26 │ 3.51837 │ └────────────┴─────────────────┴──────────────────┴───────────────────┴──────────────────┴────────────────┴──────────────────┘