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)