Introduction¶
In this notebook we solve the computational exercises of the first problem set in macroeconomics 1, due by the 16th of November. This document is structured as follows: First, we implement functions necessary to compute the analytical solution, which boils down to two functions value_function_analytical and policy_function_analytical which map capital levels to the value und policy level, respectively, given the analytical solution. Second, we implement the functions necessary for the numerical implementation. The main function in this section is value_func_iteration and _internal_value_func_iteration which provide a wrapper and implementation of the value function iteration algorithm. At last we define a few helper functions to present the solutions to exercises in a nice fashion.
If you are only interested in the solutions please go directly to subsection "Exercise Solution".
Preliminaries¶
You have to install numpy, pandas, seaborn, matplotlib and numba. (We use numba for speed but speed gains are limited since the exercise is so simple..)
In [1]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from functools import partial
from numba import njit
sns.set_style("whitegrid")
Analytical Solution¶
In [2]:
def alpha_one(alpha, beta):
"""Alpha one as computed in exercise 1.8
Args:
alpha (float): Capital exponent in production function
F(k, l) = A * k ** alpha; 0 <= alpha.
beta (float): Discount factor in household problem; 0 < beta < 1.
Returns:
out (float): Value of alpha_0.
"""
return alpha / (1 - alpha * beta)
In [3]:
def alpha_zero(beta, technology, alpha_1):
"""Alpha zero as computed in exercise 1.8
This computes alpha zero dependent on alpha_1. In principle alpha_0
can be formulated without a dependence on alpha_1 but since I am bad
in "algebra" I avoid taking the risk of making an error.
Args:
beta (float): Discount factor in household problem; 0 < beta < 1.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
alpha_one (float): The computed coefficient using function ``alpha_one``.
Returns:
out (float): Value of alpha_0.
"""
alpha_one_beta = alpha_1 * beta
out = alpha_one_beta * np.log(alpha_one_beta) + (1 + alpha_one_beta) * np.log(
technology / (1 + alpha_one_beta)
)
out /= 1 - beta
return out
In [4]:
def value_function_analytical(capital, alpha_0, alpha_1):
"""Analytical value function guess.
Args:
capital (float or np.ndarray): Capital.
alpha_0 (float): Coefficient.
alpha_1 (float): Coefficient.
Returns:
value (float or np.ndarray): Value given capital.
"""
value = alpha_0 + alpha_1 * np.log(capital)
return value
In [5]:
def policy_function_analytical(capital, alpha_1, alpha, beta, technology):
"""Analytical policy function guess.
Args:
capital (float or np.ndarray): Capital.
alpha_1 (float): Coefficient.
alpha (float): Capital exponent in production function
F(k, l) = A * k ** alpha; 0 <= alpha.
beta (float): Discount factor in household problem; 0 < beta < 1.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
Returns:
control (float or np.ndarray): Control policy.
"""
control = alpha_1 * beta * technology * (capital ** alpha)
control /= 1 + alpha_1 * beta
return control
In [6]:
@njit
def production_func(capital, technology, alpha):
"""Production function.
F(k, l) = A * k ** alpha, where k denotes capital.
Args:
capital (float or np.ndarray): Capital.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
alpha (float): Capital exponent in production function.
F(k, l) = A * k ** alpha; 0 <= alpha.
Returns:
production (float or np.ndarray): Production level.
"""
production = technology * capital ** alpha
return production
In [7]:
@njit
def utility_func(capital_today, capital_tomorrow, alpha, technology):
"""Utility function only dependent on capital.
Utility function is given as U(c) where c is the level of consumption. The level
of consumption can be substituted by c = log(f(k) - k'), where f denotes the
production function, k the level of capital today and k' the level of capital
tomorrow.
Args:
capital_today (float or np.ndarray): Capital levels today.
capital_tomorrow (float or np.ndarray): Capital levels tomorrow.
alpha (float): Capital exponent in production function.
F(k, l) = A * k ** alpha; 0 <= alpha.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
Returns:
utility (float or np.ndarray): Utility level given "consumption".
"""
utility = np.log(
production_func(capital_today, technology, alpha) - capital_tomorrow
)
return utility
In [8]:
@njit
def sup_norm_metric(arr1, arr2):
"""Supremum norm metric defined on arrays.
Computes, with n := length(arr1),
max_{i=1,...,n} |arr1_i - arr2_i|
Args:
arr1 (np.ndarray): Array of length n.
arr2 (np.ndarray): Array of length n.
Returns:
distance (float): Supremum norm.
"""
distance = np.max(np.abs(arr1 - arr2))
return distance
In [9]:
@njit
def resize(array2d, new_size):
"""Resize two dimensional array.
Args:
array2d (np.ndarray): Two dim numpy array.
new_size (tuple): New dimensions.
Returns:
new (np.ndarray): Array with new dimensions where new
entries are filled with zeroes.
"""
new = np.zeros(new_size, array2d.dtype)
n, m = array2d.shape
new[:n, :m] = array2d
return new
In [10]:
@njit
def double_array_size(array):
"""Double array size.
Args:
array (np.ndarray): Two dim. numpy array of shape (n, m).
Returns:
array (np.ndarray): Two dim. numpy array of shape (2*n, m) with the last n rows
filled with zeros.
"""
array = array.copy()
nrows, ncols = array.shape
array = resize(array, (2 * nrows, ncols))
return array
Algorithm (Value Function Iteration):¶
In [11]:
@njit
def _internal_value_func_iteration(
alpha, beta, technology, capital_grid, tolerance, initial_array_size,
):
"""Internal value function iterations.
The algorithm implemenation works as follows. We fill an array with
columns corresponding to values in the capital grid and rows
corresponding to iterations. The first guess fills the first row.
From here on we iteratively fill each row by solving the Bellman
equation. Since the previous row of the array gives the values of
the Value function of the previous iteration this amounts to solving
a simple maximum problem over a finite array. In each iteration the
supremum distance is then computed on consecutive rows. As soon as
this distance falls below the given tolerance the algorithm terminates.
Args:
alpha (float): Capital exponent in production function.
F(k, l) = A * k ** alpha; 0 <= alpha.
beta (float): Discount factor in household problem; 0 < beta < 1.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
capital_grid (np.ndarray or list-like): List of discretized capital
values of which we optimize.
tolerance (float): Positive float indicating the tolerance for the
numerical approximation.
return_info (bool): Should extra information be returned? Defaults to
True.
initial_array_size (int): Initial size of array which is filled during
the iteration process.
"""
distance = np.inf
value_array = np.zeros((initial_array_size, len(capital_grid)))
policy_array = np.empty(len(capital_grid))
n = 1
while distance > tolerance:
# initial array size is preallocated and must be increased if necessary
if n == initial_array_size:
value_array = double_array_size(value_array)
initial_array_size *= 2
# algorithm:
for j, capital in enumerate(capital_grid):
old_values = value_array[n - 1, :]
to_max = (
utility_func(capital, capital_grid, alpha, technology)
+ beta * old_values
)
argmax = np.argmax(to_max)
value_array[n, j] = to_max[argmax]
policy_array[j] = capital_grid[argmax]
distance = sup_norm_metric(value_array[n - 1, :], value_array[n, :])
n += 1
return value_array, policy_array, n
In [12]:
def value_func_iteration(
alpha,
beta,
technology,
capital_grid,
tolerance=10e-6,
return_info=False,
initial_array_size=1_000,
):
"""Value function iteration.
For an explaination see the corresponding interal function.
Args:
alpha (float): Capital exponent in production function.
F(k, l) = A * k ** alpha; 0 <= alpha.
beta (float): Discount factor in household problem; 0 < beta < 1.
technology (float): Technology parameter of production function (A) in
F(k, l) = A * k ** alpha; 0 <= A.
capital_grid (np.ndarray or list-like): List of discretized capital
values of which we optimize.
tolerance (float): Positive float indicating the tolerance for the
numerical approximation. Defaults to 10e-6.
return_info (bool): Should extra information be returned? Defaults to
True.
initial_array_size (int): Initial size of array which is filled during
the iteration process. Defaults to 1_000.
Returns:
result (pd.DataFrame): Frame containing values given capital grid.
info (dict): Extra information on the optimization; only returned if
return_info is True.
"""
capital_grid = np.array(capital_grid)
value_array, policy_array, n = _internal_value_func_iteration(
alpha, beta, technology, capital_grid, tolerance, initial_array_size
)
info = {
"number_of_iterations": n - 1,
"value_array": value_array[:n, :],
}
result = np.c_[capital_grid, value_array[n - 1, :], policy_array]
result = pd.DataFrame(result, columns=["capital", "values", "policy"])
out = (result, info) if return_info else result
return out
Computation¶
Auxiliary function for exercise presenation¶
In [13]:
def exercise_solution(
capital_grid, alpha=0.3, beta=0.6, technology=20, return_df=False
):
"""Solution to exercises 1.9-1.12.
Plot value and policy function on capital grid with numerical and
analytical solution.
Args:
return_df (bool): If True data frame with results is returned.
Defaults to false.
Returns:
df (pd.DataFrame): Results data frame in tidy format.
"""
# analytical solution
alpha_1 = alpha_one(alpha, beta)
alpha_0 = alpha_zero(beta, technology, alpha_1)
capital_grid_fine = np.linspace(2, 10, num=1000)
analytical = {
"values": value_function_analytical(capital_grid_fine, alpha_0, alpha_1),
"policy": policy_function_analytical(
capital_grid_fine, alpha_1, alpha, beta, technology
),
}
# numerical solution
df = value_func_iteration(alpha, beta, technology, capital_grid)
# create plot
fig, ax = plt.subplots(1, 2)
fig.set_size_inches(18, 8)
for i, func_type in enumerate(["values", "policy"]):
# seaborn plots
ana = sns.lineplot(
x=capital_grid_fine,
y=analytical[func_type],
ax=ax[i],
linewidth=2,
label="analytical",
)
num = sns.lineplot(
x=df["capital"], y=df[func_type], ax=ax[i], linewidth=2, alpha=0.7, label="numerical"
)
# customize axis
ax[i].set_frame_on(False)
ax[i].set_title(f"{func_type} function".capitalize(), fontsize=16, pad=-15)
ax[i].legend(
loc="lower right", frameon=False, bbox_to_anchor=(0.96, 0.14), fontsize=14
)
ax[i].xaxis.label.set_size(12)
ax[i].yaxis.label.set_size(12)
plt.show()
# return results (for debugging)
out = df if return_df else None
return out
Exercise Solutions:¶
Exercise 1.9 and 1.10¶
In [14]:
exercise_solution(capital_grid=np.arange(2, 12, 2))
Exericse 1.11 and 1.12¶
In [15]:
# capital grid 2 to 10 with 50 steps
exercise_solution(capital_grid=np.linspace(2, 10, num=25))
In [16]:
# capital grid 2 to 10 with 100 steps
exercise_solution(capital_grid=np.linspace(2, 10, num=100))
In [17]:
# capital grid 2 to 10 with 1000 steps
exercise_solution(capital_grid=np.linspace(2, 10, num=1_000))
Exercise 1.13¶
In [18]:
alpha_1 = alpha_one(alpha=0.3, beta=0.6)
policy_func = partial(policy_function_analytical, alpha_1=alpha_1, alpha=0.3, beta=0.6, technology=20)
In [19]:
capital_path = np.zeros(21)
capital_path[0] = 2
for i in range(20):
capital_path[i+1] = policy_func(capital_path[i])
In [20]:
plt.figure(figsize=(10, 6))
sns.lineplot(x=np.arange(len(capital_path)), y=capital_path, linewidth=3)
plt.title("Optimal Capital Path with $k_0 = 2$")
plt.xlabel("Time: $t=0,\dots,20$")
plt.ylabel("Capital")
plt.xticks(np.arange(0, len(capital_path), 2))
plt.box(False)
plt.show()
