7 Aplicaciones I

Cargamos lo paquetes necesarios:

using LinearAlgebra
using DelimitedFiles
using Plots

7.1 Modelo Insumo-Producto

Aplicación tomada del libro “Matemática para Economistas con Excel y Matlab” de Alicia Bernardello y coautores

Características del modelo:

  • Describe relaciones interindustriales o intersectoriales.
  • Usa dichas relaciones para predecir cambios en la demanda para predecir causados por cambios en la demanda autónoma del los productos finales.

Suponga una economía con tres sectores:

Prod/Ins I II III
I 0.2 0.6 0
II 0.2 0 0.2
III 0.4 0.2 0.5
  • Columna 1: Por cada peso producido del producto I se usan: 0.2 pesos del mismo producto, 0.2 pesos del producto II, y 0.4 pesos del producto III.
  • La demanda de cada una de las tres industrias está compuesta de la demanda derivada de las otras industrias y la demanda autónoma:

\[x_1 = 0.2x_1+0.6x_2+d_1 \] \[ x_2 = 0.2x_1+0.2x_3+d_2 \] \[ x_3 = 0.4x_1+0.2x_2+0.5x_3+d_3 \]

  • Definamos la matriz de coeficientes técnicos como:
A = [0.2 0.6 0; 0.2 0 0.2; 0.4 0.2 0.5]
3×3 Array{Float64,2}:
 0.2  0.6  0.0
 0.2  0.0  0.2
 0.4  0.2  0.5
  • Si definimos los vectores \(x=(x_1,x_2,x_3)‘\) y \(d=(d_1,d_2,d_3)’\) tenemos:

\[x = Ax + d\]

  • Entonces:

\[x=(I-A)^{-1}d\]

  • La matriz \(I-A\) es denominada matriz de Leontief.
  • La matriz \(A\) tiene todo sus autovalores en el círculo unitario y \(I-A\) es no singular. Entonces \((I-A)^{-1}\) es no negativa.
eigA, autvalA = eigen(A);
eigA
3-element Array{Float64,1}:
  0.7445080792964658 
  0.14315036673345607
 -0.18765844602992152
autvalA
3×3 Array{Float64,2}:
 -0.365644  -0.68212     0.797453
 -0.331827   0.0646304  -0.515232
 -0.869595   0.728379   -0.314014
Identity = Matrix{Int64}(I, 3, 3)
detL = det(Identity-A)
0.26

Como esperábamos:

invL = inv(Identity-A)
3×3 Array{Float64,2}:
 1.76923   1.15385  0.461538
 0.692308  1.53846  0.615385
 1.69231   1.53846  2.61538
  • Si \(d_1=100\), \(d_2=50\), y \(d_3=200\), las demandas sectoriales serían:
d = [100; 50; 200];
x = invL*d
3-element Array{Float64,1}:
 326.9230769230769
 269.2307692307692
 769.2307692307693
  • Nos interesa calcular el efecto de un cambio en la demanda autónoma de un producto sobre la demanda total (dadas las relaciones entre sectores).

\[\Delta x = (I-A)^{-1} \Delta d\]

  • Suponga ahora \(d’_1=110\), \(d’_2=100\), y \(d’_3=180\). Entonces el cambio en la producción será:
Dd = [10; 0; 0];
Dx = invL*Dd
3-element Array{Float64,1}:
 17.692307692307693
  6.923076923076923
 16.923076923076923

7.2 Equilibrio de Mercado con Dos Bienes

Aplicación tomada del libro “Matemática para Economistas con Excel y Matlab” de Alicia Bernardello y coautores

El Modelo:

  • Función de Demanda del producto 1: \[Q_{d1} = 40 - 2 P_1 + P_2\]
  • Función de Oferta del producto 1: \[Q_{o1} = -5 + 3 P_1 - P_2\]
  • Equilibrio en el mercado 1: \[Q_{d1}=Q_{o1}=Q_1\]
  • Función de Demanda del producto 2: \[Q_{d2} = 90 + P_1 - P_2\]
  • Función de Oferta del producto 2: \[Q_{o2} = -2 + 2 P_2 \]
  • Equilibrio en el mercado 2: \[Q_{d2}=Q_{o2}=Q_2\]

El modelo puede ser resuelto usando el método de reemplazo para encontrar la solución al sistema de ecuaciones. Bastante tedioso. Definamos mas bien el siguiente vector: \[x = (Q_1,Q_2,P_1,P_2)'\]

Definamos entonces: \[ Ax = b\]

La solución del sistema de ecuaciones es entonces (si \(A\) es invertible): \[x = A^{-1} \times b\]

A = [1 0 2 -1;
     1 0 -3 1;
     0 1 -1 1;
     0 1 0 -2]
4×4 Array{Int64,2}:
 1  0   2  -1
 1  0  -3   1
 0  1  -1   1
 0  1   0  -2
b = [40; -5; 90; -2]
4-element Array{Int64,1}:
 40
 -5
 90
 -2
dA = det(A)
-13.0

Solución:

x = A \ b   # alternativamente x =inv(A)*b
4-element Array{Float64,1}:
 29.76923076923076
 75.6923076923077 
 24.53846153846154
 38.84615384615385

Entonces: \(Q_1=29.7\) ,\(P_1=24.5\), \(Q_2=75.7\), \(P_2=38.8\).

7.3 Mínimos Cuadrados Ordinarios

Construiremos dos vectores \(X_1\) y \(X_2\) de \(N 1\) cada uno y con ello un vector \(y\) de acuerdo a:

\[y_i = 3 + 0.5X_{1i} + 0.9X_{2i} + \epsilon_i\]

donde \(X_{1i} N(1,4)\), \(X_{2i} N(2,4)\) y \(N(0,1)\).

A partir de estos datos usaremos la fórmula estándar de MCO para encontrar los coeficientes estimados. \[\hat{\beta} = (X'X)^{-1} X'y\]

  • Generando Datos:
N = 100;
X1 = 1.0*ones(N,1) + 2*randn(N,1);
X2 = 2.0*ones(N,1) + 2*randn(N,1);
eps = randn(N,1);
y  = 3.0*ones(N,1) + 0.5*X1 +0.9*X2 + eps;

data = [y X1 X2]
100×3 Array{Float64,2}:
 8.01438    0.210875    2.89604 
 1.80203   -3.65555     0.654704
 2.88655   -2.67982     1.60507 
 7.04636    4.11119     2.11566 
 3.557     -1.1733      0.408679
 5.89969    2.694       2.70327 
 9.59384    2.17003     3.4518  
 5.25342    1.56164     1.88716 
 6.80268   -0.0836415   4.58659 
 4.12239    0.720062    0.893211
 5.26949    1.62999     2.21036 
 0.131776  -2.93849    -0.108095
 8.86826    5.8032      2.81857 
 ⋮                              
 5.88792    3.91547     0.619451
 8.59154    4.34637     3.55927 
 2.23504    3.3423     -1.50933 
 3.38372   -2.32858     1.00137 
 4.84357    0.646629    1.15271 
 3.73546   -2.75804     3.58588 
 6.89851    0.231349    3.4753  
 7.68679    2.00034     3.46007 
 5.60082    1.8471      1.56512 
 3.77353    2.395       2.59936 
 3.92604   -1.98861     2.1217  
 4.34771   -2.94256     3.72675

Guardamos la información en un archivo “txt”.

writedlm("testols.txt", data)

Ahora cargando los datos y estimamos los parámetros vía MCO:

impdata = readdlm("testols.txt")
100×3 Array{Float64,2}:
 8.01438    0.210875    2.89604 
 1.80203   -3.65555     0.654704
 2.88655   -2.67982     1.60507 
 7.04636    4.11119     2.11566 
 3.557     -1.1733      0.408679
 5.89969    2.694       2.70327 
 9.59384    2.17003     3.4518  
 5.25342    1.56164     1.88716 
 6.80268   -0.0836415   4.58659 
 4.12239    0.720062    0.893211
 5.26949    1.62999     2.21036 
 0.131776  -2.93849    -0.108095
 8.86826    5.8032      2.81857 
 ⋮                              
 5.88792    3.91547     0.619451
 8.59154    4.34637     3.55927 
 2.23504    3.3423     -1.50933 
 3.38372   -2.32858     1.00137 
 4.84357    0.646629    1.15271 
 3.73546   -2.75804     3.58588 
 6.89851    0.231349    3.4753  
 7.68679    2.00034     3.46007 
 5.60082    1.8471      1.56512 
 3.77353    2.395       2.59936 
 3.92604   -1.98861     2.1217  
 4.34771   -2.94256     3.72675 
N = 100
y = impdata[:,1]
X1 = impdata[:,2]
X2 = impdata[:,3]

X = [ones(N,1) X1 X2]  # Concatenamos

beta = (X'*X)\(X'*y)   # alternativa beta = inv(X'*X)*(X'*y)
3-element Array{Float64,1}:
 2.949550718656315 
 0.5329621878953407
 0.9480453633565741
  • Graficando la predicción y los errores del modelo.
err = y - X*beta
yhat = X*beta
obs = collect(1:N)

plt_yhat = plot(obs,[y yhat], xlabel="Observación", ylabel="y", title = "Ajuste MCO", 
                  color=[:blue :red], label=["Observado" "Ajustado"], legend = true, 
                  linewidth = 2, shape = [:circle :diamond], grid = true)
display(plt_yhat)

svg

plt_err = plot(obs,err, xlabel="Observación", ylabel="Error de Regresión", title = "Ajuste MCO", 
                  color=[:blue], legend = false, linewidth = 2, shape = [:circle], grid = true)
display(plt_err)

svg