```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 8: Deterministic Time Series Models II

## Nonlinear First-Order Difference Equations

Recall the Solow growth model with no exogenous productivity growth (written in "per worker" terms):

\begin{align}
y_t & = Ak_t^{\alpha}\\
k_{t+1} & = i_t + (1-\delta)k_{t}\\
y_t & = c_t + i_t\\
c_t & = (1-s)y_t,
\end{align}

where $y_t$ is output per worker, $k_t$ is capital per worker, $c_t$ is consumption per worker, and $i_t$ is investment per worker. We've seen that by treating investment as an exogenous quantity, then the capital accumulation equation could be viewed as a linear first-order difference equation. But in the Solow model, investment is *not* exogenous and is equal to the savings rate times output per worker:

\begin{align}
i_t & = sy_t = sAk_t^{\alpha}
\end{align}

Therefore the equilibrium law of motion for capital per worker is:

\begin{align}
k_{t+1} & = sAk_t^{\alpha} + (1-\delta)k_{t}, \label{eqn:capital_solution}
\end{align}

which is a *nonlinear* first-order difference equation. In equilibrium capital in period $t+1$ is a *concave* (i.e., increasing at a decreasing rate) function of capital in period $t$. We can iterate on the nonlinear law of motion just like we iterated on the linear difference equation.

Let's suppose the following values for the simulation:

| $$k_0$$ | $$s$$ | $$A$$  | $$\alpha$$ | $$\delta$$  | $$T$$  |
|---------|-------|--------|------------|-------------|--------|
| 10      | 0.1   | 10     | 0.35       |  0.1        | 101    |

Where $T$ is the total number of simulation periods (i.e., $t$ will range from $0$ to $100$).


```python
# Create variables 'k0', 's', 'A', 'alpha','delta','T' to store parameter values for the simulation


# Initialize capital per worker variable 'capital' as an array of zeros with length T


# Set first value of capital equal to k0


# Iterate over t in range(T-1) to update subsequent values in the capital array


# Construct a plot of simulated capital per worker



```

It looks like capital per worker is converging toward a steady state value near $35$. You can verify that the steady state is in fact $34.55$. 

So computing nonlinear difference equations is essentially the same as computing linear differrence equations. However, establishing the stability properties of nonlinear difference equations is more involved and we'll skip that. But in case you're interested, I've included an optional discussion at the bottom of this notebook.

## Systems of Difference Equations

Most macroeconomic models are dynamic *systems* of equations. They are collections of several equations that simultaneously determine several variables. Working with them can be a little intimidating, but in principle, it's similar to working with a single dynamic equation.

### Example: Solow Growth Model Without Exogenous Growth

The basic Solow growth model is a system of four equations that determine values of four variables $k_t$, $y_t$, $c_t$, and $i_t$. In the Solow model, capital per worker is *state variable* because the the value of $k_t$ was actually determined in period $t-1$. In fact, in the Solow model without exogenous population of productivity growth, capital is the *only* state variable. That means that if you know the value of capital per worker, then you can compute all of the other quantities. The nonstate variables are called *control* variables.

To make this point clear, here is the Solow growth model *solved* in terms of $k_t$:

\begin{align}
y_t & = Ak_t^{\alpha}\\
i_t & = sAk_t^{\alpha}\\
c_t & = (1-s)Ak_t^{\alpha}\\
k_{t+1} & = sAk_t^{\alpha} + (1-\delta)k_{t}.
\end{align}

To simulate all of the endogenous variables, do these two steps:

1. Simulate an array of values for $k_t$
2. Compute the other variables using the array of $k_t$ values.

Since we've already simulated capital, we can just compute simulated values for $y_t$, $c_t$, and $i_t$ right away.


```python
# Create variables 'output', 'consumption', and 'investment' equal to the respective variables in the model


# Construct a plot of simulated output per worker, consumption per worker, investment per worker, capital per worker.


# Create legend. Use ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) to put legend outside of axis


# Add grid if you want

```

Of course, if we want to compute multiple simulations, then we will want to write a function to compute the simulated values so that we don't have to repeat ourselves. The function in the following cell returns a `DataFrame` with columns containing simulated values of capital per worker, output per worker, consumption per worker, and investment per worker.


```python
# Define a function that returns a DataFrame of simulated value from the Solow model with no exogenous growth. CELL PROVIDED
def solow(s,A,alpha,delta,T,k0):
    
    '''Function for computing a simulation of the Solow growth model without exogenous growth.
    
    The model is assmed to be in "per-worker" quantities:
    
            y[t]   = A*k[t]^alpha
            k[t+1] = i[t] + (1-delta)*k[t]
            y[t]   = c[t] + i[t]
            c[t]   = (1-s)*y[t]
    
    Args:
        s (float):      Saving rate
        A (float):      TFP
        alpha (float):  Capital share in Cobb-Dougla production function
        delta (float):  Capital depreciation rate
        T (int):        Number of periods to simulate
        k0 (float):     Initial value of capital per worker

     
     Returns:
         Pandas DataFrame
     '''
    
    # Initialize capital per worker values
    capital = np.zeros(T)

    # Set first value of capital per worker equal to k0
    capital[0] = k0

    # Iterate over t in range(T-1) to update subsequent values in the capital array
    for t in range(T-1):
        capital[t+1] = s*A*capital[t]**alpha + (1-delta)*capital[t]
    
    # Compute the values of the other variables
    output = A*capital**alpha
    consumption = (1-s)*output
    investment = s*output
    
    # Put simulated data into a DataFrame
    df = pd.DataFrame({'output':output,'consumption':consumption,'investment':investment,'capital':capital})
    
    # Return the simulated data
    return df
```

Use the function `solow()` to simulate the Solow model with the same parameters that we've been using.


```python
# Simulate the model and store results in a variable called 'df'


# Print the first five rows of df

```

Now, use the function `solow()` to simulte the model for five different initial values of capital per worker: $k_0 = 10, 20, 30, 40, 50$. Plot the trajectories for $y_t$ together.


```python
# Create variable 'initial_ks' that store the five initial values of k[t]


# Create a figure and axis


# Iterate over k0 in initial_ks and plot


# Add title


# Create legend. Use ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) to put legend outside of axis


# Add grid if you want

```

### Solow Model with Exogenous Labor

Now consider the Solow growth model with exogenous labor growth. That is, let $L_t$ denote the quantity of labor and suppose that the population (and therefore the labor supply) grows at the constant rate $n$. This means:

\begin{align}
L_{t+1} & = (1+n)L_t,
\end{align}

where the intial value of labor $L_0$ must be given. The rest of the model in aggregate terms is:

\begin{align}
Y_t & = AK_t^{\alpha}L_t^{1-\alpha}\\
K_{t+1} & = I_t + (1-\delta)K_{t}\\
Y_t & = C_t + I_t\\
C_t & = (1-s)Y_t,
\end{align}

where the intial value of capital $K_0$ is also given.

Since capital *and* labor are both determined in the previous period by their respective laws of motion, the Solow model with exogenous labor growth has *two* state variables: $K_t$ and $L_t$. We can solve the capital law of motion for capital in terms of only the variables capital and labor and so the two equations that determine how the state of the economy evolves are:

\begin{align}
K_{t+1} & = sAK_t^{\alpha}L_t^{1-\alpha} + (1-\delta)K_{t}\\
L_{t+1} & = (1+n)L_t
\end{align}

If we iterate on these two to compute simulations of $K_t$ and $L_t$, then we can compute $Y_t$, $C_t$, and $I_t$ easily.

Let's suppose the following values for a simulation:

| $$L_0$$ | $$n$$  | $$K_0$$ | $$s$$  | $$A$$  | $$\alpha$$ | $$\delta$$ | $$T$$ |
|---------|--------|---------|--------|--------|------------|------------|-------|
| 1       | 0.01   | 10      | 0.1    | 10     | 0.35       | 0.1        | 101   |

Compute simulated values for $K_t$ and $L_t$ and use those values to compute and plot simulated series for output $Y_t$ and output *per worker* $Y_t/L_t$ side-by-side in the same figure.


```python
# Create variables 'K0', 'L0', 'n', 's', 'A', 'alpha','delta','T' to store parameter values for the simulation


# Initialize capital variable 'capital' as an array of zeros with length T


# Set first value of capital equal to K0


# Initialize the labor variable 'labor' as an array of zeros with length T


# Set first value of capital equal to L0


# Iterate over t in range(T-1) to update subsequent values in the capital and labor arrays


# Create variables 'output' and 'output_pw' that store output and output per worker


# Create figure


# Construct a plot of simulated output


# Construct a plot of simulated output per worker



```

Because of exogenous labor growth, there is long-run growth in output but not output per worker. The function in the cell below simulates the Solow model with exogenous labor growth.


```python
# Define a function that returns a DataFrame of simulated value from the Solow model with exogenous labor growth. CELL PROVIDED
def solow_w_L(s,A,alpha,delta,n,T,K0,L0):
    
    '''Function for computing a simulation of the Solow growth model with exogenous labor.
    
            Y[t]   = A*K[t]^alpha*L[t]^(1-alpha)
            K[t+1] = I[t] + (1-delta)*K[t]
            Y[t]   = C[t] + I[t]
            C[t]   = (1-s)*Y[t]
            L[t+1] = (1+n)*L[t]
    
    Args:
        s (float):      Saving rate
        A (float):      TFP
        alpha (float):  Capital share in Cobb-Douglas production function
        delta (float):  Capital depreciation rate
        n (float):      Labor growth rate
        T (int):        Number of periods to simulate
        k0 (float):     Initial value of capital per worker
        L0 (float):     Initial labor
     
     Returns:
         Pandas DataFrame
     '''
    
    # Initialize capital values
    capital = np.zeros(T)

    # Set first value of capital equal to K0
    capital[0] = K0
    
    # Initialize labor values
    labor = np.zeros(T)

    # Set first value of labor equal to L0
    labor[0] = L0

    # Iterate over t in range(T-1) to update subsequent values in the capital and labor arrays
    for t in range(T-1):
        capital[t+1] = s*A*capital[t]**alpha*labor[t]**(1-alpha) + (1-delta)*capital[t]
        labor[t+1] = (1+n)*labor[t]
    
    # Compute the values of the other aggregate variables
    output = A*capital**alpha*labor**(1-alpha)
    consumption = (1-s)*output
    investment = s*output
    
    # Compute the values of the "per worker" variables
    capital_pw = capital/labor
    output_pw = output/labor
    consumption_pw = consumption/labor
    investment_pw = investment/labor
    
    # Put simulated data into a DataFrame
    df = pd.DataFrame({'output':output,
                       'consumption':consumption,
                       'investment':investment,
                       'capital':capital,
                       'output_pw':output_pw,
                       'consumption_pw':consumption_pw,
                       'investment_pw':investment_pw,
                       'capital_pw':capital_pw})
    
    # Return the simulated data
    return df
```

Use the function `solow_w_L()` to replicate the previous simulation and plots.


```python
# Replicate the previous simulation exercise using the solow_no_exo_growth() function. CELL PROVIDED

# Simulate the model
df = solow_w_L(s,A,alpha,delta,n,T,K0,L0)

# Create a figure
fig = plt.figure(figsize=(12,4))

# Construct a plot of simulated output
ax = fig.add_subplot(1,2,1)
ax.plot(df['output'],lw=3,alpha=0.75)
ax.set_title('Ouput')
ax.grid()

# Construct a plot of simulated output per worker
ax = fig.add_subplot(1,2,2)
ax.plot(df['output_pw'],lw=3,alpha=0.75)
ax.set_title('Ouput per Worker')
ax.grid()
```

## Stability of Nonlinear First-Order Difference Equations (Optional)

Unlike the linear first-order difference equation, stability (i.e., whether the process diverges to infinity in absolute value) of nonlinear first-order difference equations is less straightforward to establish. In fact, proving or disproving *global stability* -- Will the process remain finite for *any* initial condition? -- is particularly challenging. An easier task is to establish the existience of *local stability*: that the process is stable in a *neighborhood* of a given point. The idea is to use the first-order Taylor series approximation (https://en.wikipedia.org/wiki/Taylor_series) to approximate the nonlinear equation with a linear one. An example using the Solow model will make this easier to explain.

Let $k^*$ denote the steady state of capital per worker in the Solow model with no exogenous growth. Then:

\begin{align}
k^* & = \left(\frac{sA}{\delta}\right)^{\frac{1}{1-\alpha}}
\end{align}

The first-order Taylor series approximation to the capital law of motion around steady state capital per worker is:

\begin{align}
k_{t+1} & \approx \left[ sA\left(k^*\right)^{\alpha} + (1-\delta)k^*\right] + \left[\alpha sA\left(k^*\right)^{\alpha-1} + 1-\delta\right]\left(k_t - k^*\right). \label{eqn:capital_approx}
\end{align}

Equation (\ref{eqn:capital_approx}) is a linear first-order difference equation. As long as the coefficient on $k_t$, 

\begin{align}
\alpha sA\left(k^*\right)^{\alpha-1} + 1-\delta,
\end{align}

is less than 1 in absolute value, then the approximated model is stable and we say that $k_t$ is stable in a neighborhood of the steady state $k^*$. Let's compute this coefficient for the Solow model using the same parameterization we used earlier:

| $$s$$ | $$A$$  | $$\alpha$$ | $$\delta$$ |
|-------|--------|------------|------------|
| 0.1   | 10     | 0.35       |  0.1       |


```python
# Parameters


# Computer steady state


# Compute coefficient


```

So steady state capital per worker is about 35 and the coefficient on $k_t$ in the approximation is 0.935. So the Solow model as we have parameterized it is stable in the neighborhood of $k^*$. Can you find the upper bound on $\alpha$ for which the model will *not* be stable given $A = 10$, $s = 0.1$ and $\delta = 0.1$?

An obvious question is: How good is the linear approximation? It turns out that for the capital law of motion in the Solow model, the linear approximation is really good. Let's plot the exact law of motion and the approximated law of motion for $k_t \in[0,100]$.


```python
# Create array of values for k[t]


# Compute k[t+1] exactly using the law of motion for capital per worker


# Compute k[t+1] approximately using the approximation to the law of motion for capital per worker


# Create figure


# Plot exact and approximate k[t+1] against for k[t]



```

The exact law of motion and the approximated law of motion intersect at the steady state and are hardly distinguishable elsewhere. You can see the difference between them only by zooming in on the extreme values.


```python
# Create figure


# Plot exact and approximate k[t+1] against for k[t] between 0 and 10


# Plot exact and approximate k[t+1] against for k[t] between 80 and 100



```

For $k_t$ near zero, the approximated law of motion is greater than the exact one but not by much. For $k_t$ near 100, the two laws of motion are still really close to each other. The point is that we can be confident that in the Solow growth model, capital will tend to converge toward the steady state regardless of its initial value.
