Apresentação

• Começaremos esse mini curso escolhendo um modelo para a simulação.
• É muito comum, na tradição mais mainstream utilizarmos modelos que lidam com o equilíbrio geral ou parcial dos agentes.
• Mas e se nós começassemos utilizando um modelo que apesar de ter um ponto de equilíbrio, nós nunca observaremos a convergência para ele?

Modelo de Goodwin (1967) - Publicado no livro: Socialism, Capitalism and Economic Growth: Essays Presented to Maurice Dobb (Feinstein, C. H.)

• É um modelo inspirado no sistema biológico de EDOs, chamado Modelo Lotka-Volterra ou Modelo Presa-Predador.

• A ideia geral do modelo de Goodwin é tentar explicar o ciclo econômico através da relação existente no processo de Produção entre Trabalhadores e Capitalistas.

O Modelo de Goodwin

Algumas premissas do Modelo para Conveniência de Modelagem:

1. Progresso Tecnológico Constante, \(a\). \[\frac{Y}{L}=a=a_0.e^{\alpha.t},\;\;\;\alpha>0\]

2. Crescimento da Força Laboral Constante, \(N\). \[N=N_0.e^{\beta.t},\;\;\;\beta>0 \]

3. Apenas dois fatores de Produção: Trabalho e Capital, ambos homogêneos e não específicos.

4. Todas as quantidades são medidas em termos reais e líquidas.

5. Todos Salários são consumidos e todos os lucros poupados e investidos (Kalecki).

obs.: Não há, necessariamente, coincidência entre o crescimento da força de laboral e do número de empregos. Logo, não há premissa de pleno emprego!

O Modelo de Goodwin

Premissas mais empíricas e discutíveis:

6. Uma relação constante entre capital e produto (efetivo). \[k=K/Y\]
7. Um salário real que se eleva na vizinhança do pleno emprego.

Observações:

1. A premissa (5) pode ser alterada para uma relação de proporcionalidade, mas essa mudança não altera a lógica do sistema.
2. A premissa (6) pode ser suavizada, mass traria questões estruturais mais complexas para o modelo.

O Modelo de Goodwin (Origem do Sistema)

Dada a definição de \(a=Y/L\), temos que a participação da renda dos trabalhadores no produto é dada por: \[u=\frac{w.L}{Y}=\frac{w}{a}\]

Consequentemente, o Profit Share em termos residuais é dada por: \[\pi=1-w/a\]

• A Taxa de Lucro: \[r_k=\frac{\pi.Y}{K}=\frac{\left(1-\frac{w}{a}\right).Y}{K}\]

Repare que pela premissa (5), temos que \(I\equiv S\equiv\pi.Y\). Portanto, temos que: \[r_k=\frac{I}{K}=\frac{\dot{K}}{K}=\frac{\dot{Y}}{Y}\] A última parte da equação acima, obtemos sabendo que a relação \(Y/K\) é constante, logo \(\dot{Y}/Y=\dot{K}/K\).

O Modelo de Goodwin (Origem do Sistema)

Pela premissa (1), temos: \[\ln{(Y/L)}=\ln{(a_0.e^{\alpha.t})}\] \[\ln Y-\ln L=\ln a_0 + \alpha.t\] \[\frac{\dot{Y}}{Y}-\frac{\dot{L}}{L}=\alpha \] \[\frac{\dot{L}}{L}=\frac{\dot{Y}}{Y}-\alpha \]

Como já conhecemos \(\dot{Y}/Y\), temos que: \[\frac{\dot{L}}{L}=\frac{(1-w/a)}{k}-\alpha\]

O Modelo de Goodwin (Origem do Sistema)

Definindo a taxa de emprego como sendo \(v\), temos: \[v=\frac{L}{N}\] Em termos dinâmicos: \[\frac{\dot{v}}{v}=\frac{\dot{L}}{L}-\frac{\dot{N}}{N} \]

Fazendo uso da premissa (2): \(\dot{N}/N=\beta\), temos:

\[\frac{\dot{v}}{v}=\frac{(1-u)}{k}-(\alpha+\beta)\]

Reorganizando:

\[\dot{v}= \left[ \frac{1}{k}-(\alpha+\beta)-\frac{u}{k}\right].v\]

O Modelo de Goodwin (Origem do Sistema)

Agora faremos uso da premissa (7). [Curva de Phillips]

Assume-se uma relação positiva entre a taxa de crescimento dos salários em termos reais e a taxa de emprego.

Fazendo uma aproximação linear da figura acima, temos:

\[\frac{\dot{w}}{w}=-\gamma+\rho.v\]

Sabendo que \(u=w/a\), temos em termos dinâmicos:

\[\frac{\dot{u}}{u}=\frac{\dot{w}}{w}-\alpha \]

O que substituindo uma na outra, gera:

\[\frac{\dot{u}}{u}=-(\gamma+\alpha)+\rho.v \]

O Modelo de Goodwin (Sistema Bi-Dimensional)

Juntando a equação de movimento de \(v\) com a de \(u\), temos então o sistema de EDOs Bi-Dimensional:

\[\dot{v}=\left[\frac{1}{k}-(\alpha+\beta)-\frac{u}{v}\right].v\;\;\;\;\;\;\;\;\;(1)\]

\[\dot{u}=\left[-(\alpha+\gamma)+\rho.v\right].u\;\;\;\;\;\;\;\;\;(2)\]

Repare que se definirmos os parâmetros como: \(a_1=\frac{1}{k}-(\alpha+\beta)\); \(b_1 = 1/k\); \(a_2 = \alpha + \gamma\); \(b2 = \rho\);

Podemos reescrever o sistema bi-dimensional como:

\[\dot{v}=(a_1-b_1.u).v\] \[\dot{u}=-(a_2-b_2.v).u\]

Que é exatamente o sistema de presa-predador de Lotka-Volterra.

O Modelo de Goodwin - Solução Analítica

• Iremos fazer uso do último sistema para resolvê-lo analiticamente.

\[\dot{v}=(a_1-b_1.u).v\] \[\dot{u}=-(a_2-b_2.v).u\]

Parte 1

• Começamos por algumas manipulações algébricas. Primeiro, multiplicamos a primeira equação por \(a_2/v\) e a segunda equação por \(a_1/u\) em ambos os lados:

\[a_2.\frac{\dot{v}}{v}=(a_1.a_2-b_1.a_2.u)\]

\[a_1.\frac{\dot{u}}{u}=-(a_1.a_2-a_1.b_2.v)\]

• Somando as equações, temos:

\[a_2.\frac{\dot{v}}{v} + a_1.\frac{\dot{u}}{u}=-b_1.a_2.u+a_1.b_2.v \]

• Que pode ser substituída em:

\[a_2.\frac{d \ln{v}}{dt} + a_1.\frac{d \ln{u}}{dt}=-b_1.a_2.u+a_1.b_2.v \]

O Modelo de Goodwin - Solução Analítica

Parte 2

• Agora, voltamos novamente para as duas EDOs e multiplicaremos a primeira por \(b_2\) e a segunda por \(b_1\) em ambos os lados:

\[\dot{v}=(a_1-b_1.u).v\] \[\dot{u}=-(a_2-b_2.v).u\]

• Originando:

\[b_2.\dot{v}=(b_2.a_1-b_2.b_1.u).v\] \[b_1.\dot{u}=-(b_1.a_2-b_1.b_2.v).u\]

• Somando as equações, temos:

\[b_2.\dot{v}+b_1.\dot{u}=(b_2.a_1.v-b_2.b_1.u.v)-(b_1.a_2.u-b_1.b_2.v.u)\] \[b_2.\dot{v}+b_1.\dot{u}=b_2.a_1.v-b_1.a_2.u\]

O Modelo de Goodwin - Solução Analítica

• Devemos reparar que o lado direito da última equação da Parte 1 é igual ao lado direito da última equação da Parte 2.

\[a_2.\frac{d \ln{v}}{dt} + a_1.\frac{d \ln{u}}{dt}=-b_1.a_2.u+a_1.b_2.v \] \[b_2.\dot{v}+b_1.\dot{u}=b_2.a_1.v-b_1.a_2.u\]

• Portanto, necessariamente como consequência, temos que ao subtrair as duas equações:

\[b_2.\dot{v}+b_1.\dot{u}-a_2.\frac{d \ln{v}}{dt} - a_1.\frac{d \ln{u}}{dt}=0\]

• Agora integrando ambos os lados em \(dt\), temos:

\[\int b_2.\dot{v} dt+\int b_1.\dot{u}dt-\int a_2.\frac{d \ln{v}}{dt}dt - \int a_1.\frac{d \ln{u}}{dt}dt=\int0dt\]

\[b_2.v+b_1.u- a_2.\ln{v} - a_1.\ln{u}=A\]

O Modelo de Goodwin - Solução Analítica

• Com o resultado anterior, podemos exponencializar e obter:

\[b_2.v+b_1.u- a_2.\ln{v} - a_1.\ln{u}=A\]

\[e^{b_2.v+b_1.u- a_2.\ln{v} - a_1.\ln{u}}=e^A\] \[e^{b_2.v}.e^{b_1.u}.e^{- a_2.\ln{v}}.e^{ - a_1.\ln{u}}=e^A\]

\[e^{b_2.v}.e^{b_1.u}.v^{- a_2}.u^{ - a_1}=B\] \[v^{- a_2}.e^{b_2.v}=B.e^{-b_1.u}.u^{a_1}\]

• Chamando os termos de:

\[X_1 = X_1(v)=v^{- a_2}.e^{b_2.v}\] \[X_2=X_2(u)=e^{-b_1.u}.u^{a_1}\]

• Temos a seguinte relação que se esbelece entre as integrais: \[X_1=B.X_2\]

O Modelo de Goodwin - Solução Analítica

• Como nosso objetivo agora é obter as curvas integrais, iremos estudar as funçoes \(X_1(v)\) e \(X_2(u)\):

\[\frac{dX_1}{dv}=-a_2.v^{-a_2-1}.e^{-b_1.u}+b_2.v^{-a_2}.e^{b_2.v}\] \[\frac{dX_1}{dv}=v^{-a_2}.e^{b_2.v}.\left(b_2-\frac{a_2}{v}\right)\] \[\frac{dX_1}{dv}=X_1.\left(b_2-\frac{a_2}{v}\right)\]

rm(list=ls())

v = seq(0.01, 12, 0.01)
a2 = 0.1
b2 = 0.05

min = a2/b2

X1 = v^(-a2)*exp(b2*v)

plot(v,X1, 
     type = "l", 
     col = "red", 
     lwd = 2, 
     main = "X1 x v - Modelo de Goodwin")

O Modelo de Goodwin - Solução Analítica

(Continuando)

\[\frac{dX_2}{du}=a_1.u^{a_1-1}.e^{-b_1.u}-b_1.v^{a_1}.e^{b_1.u}\] \[\frac{dX_2}{du}=u^{a_1}.e^{-b_1.u}.\left(\frac{a_1}{u}-b_1\right)\] \[\frac{dX_2}{du}=X_2.\left(\frac{a_1}{u}-b_1\right)\]

u = seq(0.01, 12, 0.01)
a1 = 0.1
b1 = 0.05

X2 = u^a1 * exp(-b1*u)

plot(u,X2, type = "l", col = "red", lwd = 2, main = "X2 x u - Modelo de Goodwin")

O Modelo de Goodwin - Curvas de Nível

• Juntando as curvas e plotando em curva de nível através da função implícita.

rm(list=ls())


alpha = 0.05 # Taxa de Crescimento da Produtividade do Trabalho
beta  = 0.04 # Taxa de Crescimento Populacional
gamma = 0.05 # Taxa de Crescimento Autônomo dos Salários Reais
rho   = 0.30 # Sensibilidade do Crescimento dos Salários ao Emprego
k     = 4.00 # Relação K/Y

a2 = alpha + gamma
b2 = rho
a1 = 1/k - (alpha+beta)
b1 = 1/k

v = seq(0.01, 1.5, 0.01)
u = seq(0.01, 1.5, 0.01)

funcao <- function(v,u){v^(-a2)*exp(b2*v)-u^a1*exp(-b1*u)}
z <- outer(v,u,funcao)

contour(v,u,z, col = "red", 
        main = "Curvas de Nível - Modelo de Goodwin",
        ylab = "u - Wage Share", xlab = "v - Employment Level")

O Modelo de Goodwin - Supercíe 3D da Função Implícita

• Uma outra forma de visualizar é transformar a função implícita \(F(v,u)\) em uma superfície tridimensional.

• Entre duas curvas de nível qualquer, temos sempre outras infinitas curvas intermediárias.

• Cada curva interna, mostra uma possibilidade de órbita fechada para um par ordenado \(v,u\).

• Script no R para a superfície em três dimensões da Função Implícita:

# Pacote para Plots em 3D.
library(plotly)

fig <- plot_ly(z = z, x = v, y = u)
fig <- fig %>% 
  add_surface(opacity = .7) %>%
  layout(
    scene=list(
      xaxis=list(title="Employment Level"),
      yaxis=list(title="Wage Share"),
      zaxis=list(title="Superficie da Funcao Implicita")
    ))

fig

O Modelo de Goodwin - Simulação Computacional (Pacote PhaseR)

• Pacote PhaseR é um pacote para plotar o digrama de fase de sistemas uni ou bi-dimensionais de EDOs.

• É um pacote que tem as funcionalidades de:

  1. Apresentar o Diagrama de Fase.
  2. Plotar as Isóclinas do Sistema.
  3. Plotar as trajetórias temporais através do cálculo por integração numérica para uma dada condição inicial da EDO.

• Antes de usar as funcionalidades do pacote, devemos criar uma função que escreve o sistema de EDOs utilizando como inputs os parâmetros, o valor da variável(is) estado(s) e o próprio tempo.

rm(list=ls())


#install.packages("phaseR")
library(phaseR)

goodwin <- function(t, y, parameters){
  
  # Definição dos Parâmetros e das Variáveis Estado
  alpha <- parameters[1]
  beta  <- parameters[2]
  gamma <- parameters[3]
  rho   <- parameters[4]
  sigma <- parameters[5]
  v     <- y[1]
  u     <- y[2]
  
  # As EDOs
  dv_dt <- (1/sigma-(alpha+beta)-u/sigma)*v
  du_dt <- (-(alpha+gamma)+rho*v)*u
  
  # A lista de quem é o output da função.
  list(c(dv_dt,du_dt))
}

• Na sequência, vamos definir os valores dos parâmetros, a dimensão do eixo x (variável v), a dimensão do eixo y (variável u) e os valores iniciais para simulação das trajetórias.

parameters = c(0.05,0.04,0.05,0.3,4)
xlim <- c(0,2.5)
ylim <- c(0,2.5)
y0 <- matrix(c(1, 1, 
               0.5, 0.5, 
               2, 1, 
               1, 2), ncol = 2, nrow = 4, byrow = TRUE)

O Modelo de Goodwin - Simulação Computacional (Pacote PhaseR)

• Agora iremos utilizar as funções flowField() para plotar o diagrama de fases.

• A função nullclines() para plotar as isóclinas no diagrama de fases.

• A função trajectories() para calcular e plotar a trajetória para uma dada (ou algumas) condição(es) inicial(ais).

goodwin.flowField <- flowField(goodwin, 
                               xlim = xlim, 
                               ylim = ylim,
                               parameters = parameters, 
                               points = 31,
                               system = "two.dim", 
                               add = FALSE, 
                               xlab = "v",
                               ylab = "u",
                               main = "Trajetórias e Diagrama de Fase")
grid()

goodwin.nullclines <- nullclines(goodwin, 
                                 xlim = xlim, 
                                 ylim = ylim,
                                 parameters = parameters, 
                                 system = "two.dim", 
                                 add.legend = FALSE, 
                                 lwd =3)

goodwin.trajectory <- trajectory(goodwin, 
                                 y0 = y0, 
                                 tlim = c(0,100),
                                 tstep = 0.01,
                                 parameters = parameters, 
                                 system = "two.dim", 
                                 col = rep("red", 4), 
                                 lwd =1.5)

O Modelo de Goodwin - Simulação Computacional (Pacote deSolve)

• Um pacote mais flexível, porém sem o diagrama de fase é o deSolve.

rm(list=ls())

library(deSolve)   # pacote para resolver as EDOs
library(latex2exp) # pacote para inserir caracteres em latex no gráfico

# parte de definição do sistema de EDOs de Goodwin
derivadas <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
      v_dot = (1/sigma-(alpha+beta)-u/sigma)*v 
      u_dot = (-(alpha+gamma)+rho*v)*u
      
    return(list(c(v_dot,u_dot)))
  })
}

• Definição dos valores dos parâmetros.

parms <- list(alpha   = 0.15,
              beta    = 0.15,
              gamma   = 0.05,
              rho     = 0.30,
              sigma   = 1)

attach(parms)
a1 = 1/sigma - (alpha+beta)
b1 = 1/sigma
a2 = (alpha + gamma)
b2 = rho

# Cálculo do centro de gravitação.
centro_1 = a1/b1
centro_2 = a2/b2

detach(parms)

# valor inicial das EDOS.
state <- c(v = 0.60, u = 0.63)

O Modelo de Goodwin - Simulação Computacional (Pacote deSolve)

• Resolução do sistema através de integração numérica pelo método de Runge-Kutta de 4ª Ordem.

• Plot das figuras em fase e contra o tempo. Além disso, o cálculo do sentido no plano de fase.

Out <- ode(y = state, times=seq(0, 50, 0.25), derivadas, parms, 
           method = "rk4")

plot(Out[,1],Out[,2], type = "l", xlab = "Tempo", main = "Wage Share and Employment Rate - Goodwin Model",
     ylab = "Wage Share and Employment Rate", col = "darkcyan", lwd = 2)
lines(Out[,1],Out[,3], type = "l", col = "tomato", lwd = 2)

plot(Out[,2],Out[,3], type = "l", 
     xlab = "Worker's share", 
     ylab = "Employment level",
     main = "Worker's Share x Employment Level - Goodwin Model", lwd =2)
points(centro_2,centro_1, lwd = 2, col = "darkviolet")
abline(h = centro_1, lty = 2, col = "grey")
abline(v = centro_2, lty = 2, col = "grey")

s <- seq(length(Out[,1] )- 1)
arrows(Out[s,2], Out[s,3], Out[s+1,2], Out[s+1,3], col = "red")

Choques no Modelo

• Umas das maneiras de estudar as propriedades do modelo é aplicando choques.
• Iremos aplicar dois. Um sobre as variáveis estado e outro sobre o parâmetro da curva de Phillips.

Choque 1 \(\uparrow (v,u)\) \((v+0.05,u+0.05)\):

parms <- list(alpha   = 0.15,
              beta    = 0.15,
              gamma   = 0.05,
              rho     = 0.30,
              sigma   = 1)

attach(parms)
a1 = 1/sigma - (alpha+beta)
b1 = 1/sigma
a2 = (alpha + gamma)
b2 = rho

centro_11 = a1/b1
centro_21 = a2/b2
Period_1 <- 2*pi/(((alpha+gamma)*(1/sigma - (alpha+beta)))^0.5)

detach(parms)

state <- c(v = 0.60, u = 0.63)

Out <- ode(y = state, times=seq(0, 50, 0.25), derivadas, parms, 
           method = "rk4")

state2 <- c(Out[nrow(Out),2]+0.05,Out[nrow(Out),3]+0.05)

parms2 <- parms

attach(parms2)
a1 = 1/sigma - (alpha+beta)
b1 = 1/sigma
a2 = (alpha + gamma)
b2 = rho

centro_12 = a1/b1
centro_22 = a2/b2

detach(parms2)


Out2 <- ode(y = state2, times=seq(50, 100, 0.25), derivadas, parms2, 
           method = "rk4")

Out_Final <- rbind(Out,Out2)


plot(Out_Final[,1],Out_Final[,2], type = "l", xlab = "Tempo", main = "Wage Share and Employment Rate - Goodwin Model",
     ylab = "Wage Share and Employment Rate", col = "darkcyan", lwd = 2)
lines(Out_Final[,1],Out_Final[,3], type = "l", col = "tomato", lwd = 2)

plot(Out_Final[,2],Out_Final[,3], type = "l", 
     xlab = "Worker's share", 
     ylab = "Employment level",
     main = "Worker's Share x Employment Level - Goodwin Model", lwd =2)
points(centro_21,centro_11, lwd = 2, col = "darkviolet")
points(centro_22,centro_12, lwd = 2, col = "darkviolet")

abline(h = centro_11, lty = 2, col = "grey")
abline(v = centro_21, lty = 2, col = "grey")

abline(h = centro_12, lty = 2, col = "grey")
abline(v = centro_22, lty = 2, col = "grey")

s <- seq(length(Out_Final[,1] )- 1)
arrows(Out_Final[s,2], Out_Final[s,3], Out_Final[s+1,2], Out_Final[s+1,3], col = "red")

Choques no Modelo

Choque 2 \(\downarrow \rho\) \(0.30\to 0.25\):

parms <- list(alpha   = 0.15,
              beta    = 0.15,
              gamma   = 0.05,
              rho     = 0.30,
              sigma   = 1)

attach(parms)
a1 = 1/sigma - (alpha+beta)
b1 = 1/sigma
a2 = (alpha + gamma)
b2 = rho

centro_11 = a1/b1
centro_21 = a2/b2
Period_1 <- 2*pi/(((alpha+gamma)*(1/sigma - (alpha+beta)))^0.5)

detach(parms)

state <- c(v = 0.60, u = 0.63)

Out <- ode(y = state, times=seq(0, 50, 0.25), derivadas, parms, 
           method = "rk4")

parms2     <- parms
parms2$rho <- 0.25

attach(parms2)
a1 = 1/sigma - (alpha+beta)
b1 = 1/sigma
a2 = (alpha + gamma)
b2 = rho

centro_12 = a1/b1
centro_22 = a2/b2

detach(parms2)


Out2 <- ode(y = state2, times=seq(50, 100, 0.25), derivadas, parms2, 
           method = "rk4")

Out_Final <- rbind(Out,Out2)


plot(Out_Final[,1],Out_Final[,2], type = "l", xlab = "Tempo", main = "Wage Share and Employment Rate - Goodwin Model",
     ylab = "Wage Share and Employment Rate", col = "darkcyan", lwd = 2)
lines(Out_Final[,1],Out_Final[,3], type = "l", col = "tomato", lwd = 2)

plot(Out_Final[,2],Out_Final[,3], type = "l", 
     xlab = "Worker's share", 
     ylab = "Employment level",
     main = "Worker's Share x Employment Level - Goodwin Model", lwd =2)
points(centro_21,centro_11, lwd = 2, col = "darkviolet")
points(centro_22,centro_12, lwd = 2, col = "darkviolet")

abline(h = centro_11, lty = 2, col = "grey")
abline(v = centro_21, lty = 2, col = "grey")

abline(h = centro_12, lty = 2, col = "grey")
abline(v = centro_22, lty = 2, col = "grey")

s <- seq(length(Out_Final[,1] )- 1)
arrows(Out_Final[s,2], Out_Final[s,3], Out_Final[s+1,2], Out_Final[s+1,3], col = "red")