Modelo SIM - Godley e Lavoie (2007) - Notação Matricial

Introdução

Em post anterior, verificamos como seria possível simular o modelo SFC do capítulo 1 do livro do Godley e Lavoie (2007). O intuito do post de hoje é mostrar que o mesmo modelo pode ser trabalhado em uma versão matricial. Dessa forma, reescreveremos o sistema de equações de modo que seja possível ter a seguinte estrutura: \[\mathbf{A.X_{t}=B.X_{t-1}+C}\] onde \(\mathbf{A}\) é a matriz \(n\times n\) de coefientes que relacionam as variáveis contemporâneas; \(\mathbf{X_{t}}\) é o vetor \(n\times 1\) das variáveis contemporâneas; \(\mathbf{B}\) é a matriz \(n\times n\) de coeficientes que relacionam as variáveis defasadas de 1ª ordem; \(\mathbf{X_{t-1}}\) é o vetor \(n\times 1\) das variáveis defasadas de 1ª ordem; \(\mathbf{C}\) é o vetor \(n\times 1\) de variáveis exógenas do modelo.

Em suma, estamos transformando o sistema de equações em um problema matricial de equações à diferenças de 1ª ordem não homogêneo.

Resolução Matricial

O sistema acima pode ser reescrito da seguinte forma:

\[\mathbf{\underset{I}{\underbrace{A^{-1}.A}}.X_{t}=A^{-1}.B.X_{t-1}+A^{-1}.C}\] O que simplificando, temos: \[\mathbf{X_t=\phi_1.X_{t-1}+\phi_2}\]

onde \(\mathbf{\phi_1=A^{-1}.B}\), \(\mathbf{\phi_2=A^{-1}.C}\) e por definição \(\mathbf{I=A^{-1}.A}\).

Podemos encontrar os pontos fixos (ou valores de estado estacionário) assumindo que \(\mathbf{X_t=X_{t-1}=X^*}\). Dessa maneira, temos que: \[\mathbf{X^*=(I-\phi_1)^{-1}.\phi_2}\] Com relação a estabilidade do sistema de equações em diferenças, necessitaremos ao fim que todos os autovalores (parte real ou imaginária) da matriz de transição \(\mathbf{A}\) possuam valor absoluto menor do que 1.

Isso posto, reorganizaremos as equações do modelo SIM para depois inputá-las na forma matricial.

Equações do Modelo SIM reescritas

As 11 equações do modelo reescritas organizando os termos ficam da seguinte forma:

\[C_s-C_d=0    (1)\] \[G_s=\overline{G}_{d}     (2)\] \[T_s-T_d=0    (3)\] \[N_s-N_d=0    (4)\] \[Y_d-W.N_s+T_s=0 (5)\] \[T_d-\theta.W.N_s=0   (6)\] \[C_d-\alpha_1.Y_d=\alpha_2.Hh_{-1} (7)\] \[Hs+T_d=Hs_{-1}+\overline{G}_{d} (8)\] \[Hh-Y_d+C_d=Hh_{-1} (9)\] \[Y-C_s-G_s=0   (10)\] \[N_d-(1/W).Y=0  (11)\]

A correspondência das variáveis se encontra no post anterior e aqui apenas lembro que a única variável exógena é \(\overline{G}_{d}\). Além disso, há uma equação escondida que será fundamental para a resoluçao do modelo, sendo ela \(Hh=Hs\). É simples e trivial porém necessária. Resgataremos isso lá na frente.

O Modelo na forma Matricial

Com a finalidade didática, abriremos abaixo as matrizes e vetores que representam nosso sistema de equações à diferença.

Em notação matricial, temos: \(\mathbf{A.X_{t}=B.X_{t-1}+C}\)

Inputando os elementos das equações, temos:

\[\begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & -1\\ 0 & 0 & 1 & -W & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -\theta.W & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -\alpha_{1} & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0 & 1 & 0 & 0\\ -1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1/W & 1 \end{bmatrix}.\begin{bmatrix}C_{s}\\ G_{s}\\ T_{s}\\ N_{s}\\ Y_{d}\\ T_{d}\\ C_{d}\\ Hs\\ Hh\\ Y\\ N_{d} \end{bmatrix}=\] \[\begin{bmatrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \alpha_{2} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}.\begin{bmatrix}C_{s-1}\\ G_{s-1}\\ T_{s-1}\\ N_{s-1}\\ Y_{d-1}\\ T_{d-1}\\ C_{d-1}\\ Hs_{-1}\\ Hh_{-1}\\ Y_{-1}\\ N_{d-1} \end{bmatrix}+\begin{bmatrix}0\\ \overline{G}_{d}\\ 0\\ 0\\ 0\\ 0\\ 0\\ \overline{G}_{d}\\ 0\\ 0\\ 0 \end{bmatrix}\]

E aqui cabe uma informação muito importante. Como a matriz \(\mathbf{A}\) será invertida durante o processo de resolução, necessitamos que ela não seja singular em um primeiro momento e nem resulte em uma matriz singular ao ser invertida. Isso porque também necessitaremos calcular \(\mathbf{(I-\phi_1)^{-1}}\).

Aqui surge então o porém mencionado anteriormente, que envolveria fazer uso da equação escondida \(Hh=Hs\).

Se inputássemos na matriz \(\mathbf{B}\), o elemento \(\mathbf{B}[8,8]=1\) ao invés de \(\mathbf{B}[8,9]=1\), obteríamos no cálculo da inversa uma matriz singular. Em termos da notação das equações, se escrevermos que \(Hs+T_d=Hs_{-1}+\overline{G}_{d}\), o sistema não pode ser solucionado. Para ser, precisamos escrever que \(Hs+T_d=Hh_{-1}+\overline{G}_{d}\), fazendo uso da equação omitida.

Resolução em Linguagem R

Apresentado o problema e a conversão de notação, vamos agora para o script que nos fornecerá os valores de estado estacionário, bem como se o sistema é estável ou não para um dado conjunto de parâmetros.

Os pacotes utilizados são o matlib e o kableExtra. O primeiro pacote é para trabalhar com cálculo matricial e o segundo para gerar tabelas em HTML. Começamos limpando o Global Environment com o comando rm(list=ls()). Na sequência, definimos os valores dos parâmetros e variáveis exógenas.

rm(list=ls())

library(matlib)
library(kableExtra)

alpha_1 <- 0.6
alpha_2 <- 0.4
W       <- 1.0
theta   <- 0.2
G_d     <- 20

Precisamos agora definir a matriz \(\mathbf{A}\). Para isso, sabemos que cada elemento da diagonal principal terá valor igual a 1. Isso porque cada equação terá sua variável endógena ordenada em igual sequência no vetor \(\mathbf{X_{t}}\) (pode-se verificar isso na última seção!). Dessa forma, a estrutura geral da matriz pode ser definida por uma matriz \(\mathbf{I}\) de tamanho \(11\times 11\). Uma forma de definí-la no R é através da função diag. Os demais elementos podem ser inputados verificando as equações e alterando individualmente onde ele se encontra na matriz.

No caso da matriz \(\mathbf{B}\), a diagonal principal não possui, a princípio, valores unitários. Portanto, dessa forma o que podemos fazer é definí-la como uma matriz de elementos \(0\) de tamanho \(11\times 11\) e depois inputarmos os elementos distintos nas respectivas linhas e colunas. O que é feito também vale para a matriz \(\mathbf{C}\), sendo essa uma matriz \(11 \times 1\). Por fim, criamos também uma matriz identidade que nos ajudará no final.

# Matriz A
A <- diag(11)

A[1,7]  <- -1
A[3,6]  <- -1
A[4,11] <- -1
A[5,4]  <- -W
A[5,3]  <- 1
A[6,4]  <- -theta*W
A[7,5]  <- -alpha_1
A[8,6]  <- 1
A[9,5]  <- -1
A[9,7]  <- 1
A[10,1] <- -1
A[10,2] <- -1
A[11,10] <- -1/W


## Matriz B
B <- matrix(0, 11, 11)

### Inputando os Elementos de B
B[7,9] <- alpha_2
B[8,9] <- 1
B[9,9] <- 1

# Observação importante. Tive que colocar B[8,9] como igual a 1.
# Não havia nenhuma equação que apontasse para Hh = Hs. Isso tornava
# a matriz singular e por consequência, não conseguíamos calcular a
# matriz inversa.

## Matriz C
C <- matrix(0,11)

C[2,1] <- G_d
C[8,1] <- G_d

I <- diag(11)

Na parte seguinte, criamos as matrizes \(\mathbf{\phi_1}\) através da multiplicação de \(\mathbf{A^{-1}.B}\) e \(\mathbf{\phi_2}\) através da multiplicação de \(\mathbf{A^{-1}.C}\). Sobre esse último termo há uma nota importante a ser feita. O resultado dela representa o vetor que fornece o equilíbrio de curto prazo das variáveis endógenas do modelo.

Por fim, temos a resolução do sistema \(\mathbf{X^*=(I-\phi_1)^{-1}.\phi_2}\) que nos fornece os vetor das variáveis no estado estacionário. Ele é apresentado em uma tabela através da função kable. Além disso, pedimos o cálculo dos autovalores da matriz de transição \(\mathbf{A}\). Como pode-se notar, entre os 11 valores possíveis, temos apenas 1 (porque o modelo é de equilíbrio único) e com valor absoluto \(<1\). Isso garante que o sistema é estável. Para que não estável fosse, o sistema teria que ter valor absoluto maior ou igual 1.

phi_1 <- inv(A) %*% B  
phi_2 <- inv(A) %*% C  

# Equilíbrio de Estado Estacionário
X_ast <- inv(I-phi_1) %*% phi_2
X_ast <- data.frame(X_ast)
rownames(X_ast) <- c("Cs*","Gs*","Ts*","Ns*","Yd*","Td*","Cd*","Hs*","Hh*","Y*","Nd*")
colnames(X_ast) <- "Valor de Estado Estacionário"

kable(X_ast) %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

DF_AV <- data.frame(eigen(phi_1)$values)
rownames(DF_AV) <- c("Autovalor 1", "Autovalor 2", "Autovalor 3", "Autovalor 4", 
                     "Autovalor 5", "Autovalor 6", "Autovalor 7", "Autovalor 8",
                     "Autovalor 9", "Autovalor 10", "Autovalor 11")
colnames(DF_AV) <- "Valor Obtido"

kable(DF_AV) %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

A título de exercício, o leitor pode verificar o que ocorre com a estabilidade do modelo quando \(\alpha_2=0\)

Referências

DE LA FUENTE, A. Mathematical Methods and Models for Economists. Cambridge University Press. 2000.

GANDOLFO, G. Economic Dynamics. Fourth Edition. Springer. Berlin. 2009.

GODLEY, W.; LAVOIE, M. Monetary Economics. An Integrated Approach to Credit, Money, Income, Production and Wealth. Palgrave Macmillan. Nova Iorque, EUA. 2007.