Quantecon Python Advanced
Quantecon Python Advanced
Python
II LQ Control 75
5 Information and Consumption Smoothing 77
i
5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 Two Representations of One Nonfinancial Income Process . . . . . . . . . . . . . . . . . . . . . . . 78
5.3 Application of Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4 News Shocks and Less Informative Shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.5 Representation of 𝜖𝑡 Shock in Terms of Future 𝑦𝑡 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.6 Representation in Terms of 𝑎𝑡 Shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.7 Permanent Income Consumption-Smoothing Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.8 State Space Representations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.9 Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.10 Simulating Income Process and Two Associated Shock Processes . . . . . . . . . . . . . . . . . . . . 89
5.11 Calculating Innovations in Another Way . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.12 Another Invertibility Issue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
ii
11.4 Better Representation of Roll-Over Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222
iii
17.8 Cattle Cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382
17.9 Models of Occupational Choice and Pay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383
17.10 Permanent Income Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384
17.11 Gorman Heterogeneous Households . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385
17.12 Non-Gorman Heterogeneous Households . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 386
iv
26 Etymology of Entropy 469
26.1 Information Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 469
26.2 A Measure of Unpredictability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 470
26.3 Mathematical Properties of Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 470
26.4 Conditional Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 471
26.5 Independence as Maximum Conditional Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 471
26.6 Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472
26.7 Statistical Divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472
26.8 Continuous distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472
26.9 Relative entropy and Gaussian distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
26.10 Von Neumann Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
26.11 Backus-Chernov-Zin Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474
26.12 Wiener-Kolmogorov Prediction Error Formula as Entropy . . . . . . . . . . . . . . . . . . . . . . . . 474
26.13 Multivariate Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475
26.14 Frequency Domain Robust Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475
26.15 Relative Entropy for a Continuous Random Variable . . . . . . . . . . . . . . . . . . . . . . . . . . . 476
27 Robustness 479
27.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 479
27.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 482
27.3 Constructing More Robust Policies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483
27.4 Robustness as Outcome of a Two-Person Zero-Sum Game . . . . . . . . . . . . . . . . . . . . . . . 485
27.5 The Stochastic Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 489
27.6 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 491
27.7 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 491
27.8 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496
v
32.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583
32.2 A Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584
32.3 Finite Horizon Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585
32.4 Infinite Horizon Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 589
32.5 Undiscounted Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 591
32.6 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593
32.7 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 601
vi
37.6 Adding Views . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 678
37.7 Bayesian Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 680
37.8 Curve Decolletage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 681
37.9 Black-Litterman Recommendation as Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . 685
37.10 A Robust Control Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 687
37.11 A Robust Mean-Variance Portfolio Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 688
37.12 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 689
37.13 Special Case – IID Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 690
37.14 Dependence and Sampling Frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 690
37.15 Frequency and the Mean Estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 692
vii
42.6 A Gradient Descent Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 800
42.7 A More Structured ML Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804
42.8 Continuation Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 811
42.9 Adding Some Human Intelligence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814
42.10 What has Machine Learning Taught Us? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 820
viii
48.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 951
48.2 The Economy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 952
48.3 Long Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 953
48.4 Asymptotic Mean and Rate of Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 978
IX Other 1045
51 Troubleshooting 1047
51.1 Fixing Your Local Environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1047
51.2 Reporting an Issue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1048
52 References 1049
Bibliography 1053
Index 1061
ix
x
Advanced Quantitative Economics with Python
CONTENTS 1
Advanced Quantitative Economics with Python
– Estimation of Spectra
– Additive and Multiplicative Functionals
– Classical Control with Linear Algebra
– Classical Prediction and Filtering With Linear Algebra
– Knowing the Forecasts of Others
• Asset Pricing and Finance
– Asset Pricing II: The Lucas Asset Pricing Model
– Elementary Asset Pricing Theory
– Two Modifications of Mean-Variance Portfolio Theory
– Irrelevance of Capital Structures with Complete Markets
– Equilibrium Capital Structures with Incomplete Markets
• Dynamic Programming Squared
– Optimal Unemployment Insurance
– Stackelberg Plans
– Machine Learning a Ramsey Plan
– Time Inconsistency of Ramsey Plans
– Sustainable Plans for a Calvo Model
– Optimal Taxation with State-Contingent Debt
– Optimal Taxation without State-Contingent Debt
– Fluctuating Interest Rates Deliver Fiscal Insurance
– Fiscal Risk and Government Debt
– Competitive Equilibria of a Model of Chang
– Credible Government Policies in a Model of Chang
• Other
– Troubleshooting
– References
– Execution Statistics
2 CONTENTS
Part I
3
CHAPTER
ONE
1.1 Overview
Orthogonal projection is a cornerstone of vector space methods, with many diverse applications.
These include
• Least squares projection, also known as linear regression
• Conditional expectations for multivariate normal (Gaussian) distributions
• Gram–Schmidt orthogonalization
• QR decomposition
• Orthogonal polynomials
• etc
In this lecture, we focus on
• key ideas
• least squares regression
We’ll require the following imports:
import numpy as np
from scipy.linalg import qr
For background and foundational concepts, see our lecture on linear algebra.
For more proofs and greater theoretical detail, see A Primer in Econometric Theory.
For a complete set of proofs in a general setting, see, for example, [Roman, 2005].
For an advanced treatment of projection in the context of least squares prediction, see this book chapter.
5
Advanced Quantitative Economics with Python
Assume 𝑥, 𝑧 ∈ ℝ𝑛 .
Define ⟨𝑥, 𝑧⟩ = ∑𝑖 𝑥𝑖 𝑧𝑖 .
Recall ‖𝑥‖2 = ⟨𝑥, 𝑥⟩.
The law of cosines states that ⟨𝑥, 𝑧⟩ = ‖𝑥‖‖𝑧‖ cos(𝜃) where 𝜃 is the angle between the vectors 𝑥 and 𝑧.
When ⟨𝑥, 𝑧⟩ = 0, then cos(𝜃) = 0 and 𝑥 and 𝑧 are said to be orthogonal and we write 𝑥 ⟂ 𝑧.
𝑦 ̂ ∶= arg min ‖𝑦 − 𝑧‖
𝑧∈𝑆
Hence ‖𝑦 − 𝑧‖ ≥ ‖𝑦 − 𝑦‖,
̂ which completes the proof.
For a linear space 𝑌 and a fixed linear subspace 𝑆, we have a functional relationship
Orthogonal Complement
Let 𝑆 ⊂ ℝ𝑛 .
The orthogonal complement of 𝑆 is the linear subspace 𝑆 ⟂ that satisfies 𝑥1 ⟂ 𝑥2 for every 𝑥1 ∈ 𝑆 and 𝑥2 ∈ 𝑆 ⟂ .
Let 𝑌 be a linear space with linear subspace 𝑆 and its orthogonal complement 𝑆 ⟂ .
We write
𝑌 = 𝑆 ⊕ 𝑆⟂
to indicate that for every 𝑦 ∈ 𝑌 there is unique 𝑥1 ∈ 𝑆 and a unique 𝑥2 ∈ 𝑆 ⟂ such that 𝑦 = 𝑥1 + 𝑥2 .
Moreover, 𝑥1 = 𝐸𝑆̂ 𝑦 and 𝑥2 = 𝑦 − 𝐸𝑆̂ 𝑦.
This amounts to another version of the OPT:
Theorem. If 𝑆 is a linear subspace of ℝ𝑛 , 𝐸𝑆̂ 𝑦 = 𝑃 𝑦 and 𝐸𝑆̂ ⟂ 𝑦 = 𝑀 𝑦, then
To see this, observe that since 𝑥 ∈ span{𝑢1 , … , 𝑢𝑘 }, we can find scalars 𝛼1 , … , 𝛼𝑘 that verify
𝑘
𝑥 = ∑ 𝛼𝑗 𝑢𝑗 (1.1)
𝑗=1
When a subspace onto which we project is orthonormal, computing the projection simplifies:
Theorem If {𝑢1 , … , 𝑢𝑘 } is an orthonormal basis for 𝑆, then
𝑘
𝑃 𝑦 = ∑⟨𝑦, 𝑢𝑖 ⟩𝑢𝑖 , ∀ 𝑦 ∈ ℝ𝑛 (1.2)
𝑖=1
𝐸𝑆̂ 𝑦 = 𝑃 𝑦
𝑃 = 𝑋(𝑋 ′ 𝑋)−1 𝑋 ′
An expression of the form 𝑋𝑎 is precisely a linear combination of the columns of 𝑋 and hence an element of 𝑆.
Claim 2 is equivalent to the statement
It is common in applications to start with 𝑛 × 𝑘 matrix 𝑋 with linearly independent columns and let
𝑃 𝑦 = 𝑈 (𝑈 ′ 𝑈 )−1 𝑈 ′ 𝑦
We have recovered our earlier result about projecting onto the span of an orthonormal basis.
𝛽 ̂ ∶= (𝑋 ′ 𝑋)−1 𝑋 ′ 𝑦
𝑋 𝛽 ̂ = 𝑋(𝑋 ′ 𝑋)−1 𝑋 ′ 𝑦 = 𝑃 𝑦
Because 𝑋𝑏 ∈ span(𝑋)
If probabilities and hence 𝔼 are unknown, we cannot solve this problem directly.
However, if a sample is available, we can estimate the risk with the empirical risk:
1 𝑁
min ∑(𝑦𝑛 − 𝑓(𝑥𝑛 ))2
𝑓∈ℱ 𝑁
𝑛=1
1.6.2 Solution
𝛽 ̂ ∶= (𝑋 ′ 𝑋)−1 𝑋 ′ 𝑦
𝑦 ̂ ∶= 𝑋 𝛽 ̂ = 𝑃 𝑦
𝑢̂ ∶= 𝑦 − 𝑦 ̂ = 𝑦 − 𝑃 𝑦 = 𝑀 𝑦
Let’s return to the connection between linear independence and orthogonality touched on above.
A result of much interest is a famous algorithm for constructing orthonormal sets from linearly independent sets.
The next section gives details.
Theorem For each linearly independent set {𝑥1 , … , 𝑥𝑘 } ⊂ ℝ𝑛 , there exists an orthonormal set {𝑢1 , … , 𝑢𝑘 } with
1.7.2 QR Decomposition
The following result uses the preceding algorithm to produce a useful decomposition.
Theorem If 𝑋 is 𝑛 × 𝑘 with linearly independent columns, then there exists a factorization 𝑋 = 𝑄𝑅 where
• 𝑅 is 𝑘 × 𝑘, upper triangular, and nonsingular
• 𝑄 is 𝑛 × 𝑘 with orthonormal columns
Proof sketch: Let
• 𝑥𝑗 ∶= col𝑗 (𝑋)
• {𝑢1 , … , 𝑢𝑘 } be orthonormal with the same span as {𝑥1 , … , 𝑥𝑘 } (to be constructed using Gram–Schmidt)
• 𝑄 be formed from cols 𝑢𝑖
Since 𝑥𝑗 ∈ span{𝑢1 , … , 𝑢𝑗 }, we have
𝑗
𝑥𝑗 = ∑⟨𝑢𝑖 , 𝑥𝑗 ⟩𝑢𝑖 for 𝑗 = 1, … , 𝑘
𝑖=1
For matrices 𝑋 and 𝑦 that overdetermine 𝛽 in the linear equation system 𝑦 = 𝑋𝛽, we found the least squares approximator
𝛽 ̂ = (𝑋 ′ 𝑋)−1 𝑋 ′ 𝑦.
Using the QR decomposition 𝑋 = 𝑄𝑅 gives
𝛽 ̂ = (𝑅′ 𝑄′ 𝑄𝑅)−1 𝑅′ 𝑄′ 𝑦
= (𝑅′ 𝑅)−1 𝑅′ 𝑄′ 𝑦
= 𝑅−1 (𝑅′ )−1 𝑅′ 𝑄′ 𝑦 = 𝑅−1 𝑄′ 𝑦
Numerical routines would in this case use the alternative form 𝑅𝛽 ̂ = 𝑄′ 𝑦 and back substitution.
1.8 Exercises
Exercise 1.8.1
Show that, for any linear subspace 𝑆 ⊂ ℝ𝑛 , 𝑆 ∩ 𝑆 ⟂ = {0}.
Exercise 1.8.2
Let 𝑃 = 𝑋(𝑋 ′ 𝑋)−1 𝑋 ′ and let 𝑀 = 𝐼 − 𝑃 . Show that 𝑃 and 𝑀 are both idempotent and symmetric. Can you give
any intuition as to why they should be idempotent?
Exercise 1.8.3
Using Gram-Schmidt orthogonalization, produce a linear projection of 𝑦 onto the column space of 𝑋 and verify this
using the projection matrix 𝑃 ∶= 𝑋(𝑋 ′ 𝑋)−1 𝑋 ′ and also using QR decomposition, where:
1
𝑦 ∶= ⎛
⎜ 3 ⎞⎟,
⎝ −3 ⎠
and
1 0
𝑋 ∶= ⎛
⎜ 0 −6 ⎞
⎟
⎝ 2 2 ⎠
def gram_schmidt(X):
"""
Implements Gram-Schmidt orthogonalization.
Parameters
----------
X : an n x k array with linearly independent columns
Returns
-------
U : an n x k array with orthonormal columns
"""
# Set up
n, k = X.shape
U = np.empty((n, k))
I = np.eye(n)
1.8. Exercises 19
Advanced Quantitative Economics with Python
# Normalize
U[:, i] = u / np.sqrt(np.sum(u * u))
return U
y = [1, 3, -3]
X = [[1, 0],
[0, -6],
[2, 2]]
First, let’s try projection of 𝑦 onto the column space of 𝑋 using the ordinary matrix expression:
Now let’s do the same using an orthonormal basis created from our gram_schmidt function
U = gram_schmidt(X)
U
Py2 = U @ U.T @ y
Py2
This is the same answer. So far so good. Finally, let’s try the same thing but with the basis obtained via QR decomposition:
Q, R = qr(X, mode='economic')
Q
array([[-0.4472136 , -0.13187609],
[-0. , -0.98907071],
[-0.89442719, 0.06593805]])
Py3 = Q @ Q.T @ y
Py3
1.8. Exercises 21
Advanced Quantitative Economics with Python
TWO
In addition to what’s in Anaconda, this lecture will need the following libraries:
2.1 Overview
In a previous lecture, we learned about finite Markov chains, a relatively elementary class of stochastic dynamic models.
The present lecture extends this analysis to continuous (i.e., uncountable) state Markov chains.
Most stochastic dynamic models studied by economists either fit directly into this class or can be represented as continuous
state Markov chains after minor modifications.
In this lecture, our focus will be on continuous Markov models that
• evolve in discrete-time
• are often nonlinear
The fact that we accommodate nonlinear models here is significant, because linear stochastic models have their own highly
developed toolset, as we’ll see later on.
The question that interests us most is: Given a particular stochastic dynamic model, how will the state of the system
evolve over time?
In particular,
• What happens to the distribution of the state variables?
• Is there anything we can say about the “average behavior” of these variables?
• Is there a notion of “steady state” or “long-run equilibrium” that’s applicable to the model?
– If so, how can we compute it?
Answering these questions will lead us to revisit many of the topics that occupied us in the finite state case, such as
simulation, distribution dynamics, stability, ergodicity, etc.
Note: For some people, the term “Markov chain” always refers to a process with a finite or discrete state space. We
follow the mainstream mathematical literature (e.g., [Meyn and Tweedie, 2009]) in using the term to refer to any discrete
time Markov process.
23
Advanced Quantitative Economics with Python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import lognorm, beta
from quantecon import LAE
from scipy.stats import norm, gaussian_kde
You are probably aware that some distributions can be represented by densities and some cannot.
(For example, distributions on the real numbers ℝ that put positive probability on individual points have no density
representation)
We are going to start our analysis by looking at Markov chains where the one-step transition probabilities have density
representations.
The benefit is that the density case offers a very direct parallel to the finite case in terms of notation and intuition.
Once we’ve built some intuition we’ll cover the general case.
In our lecture on finite Markov chains, we studied discrete-time Markov chains that evolve on a finite state space 𝑆.
In this setting, the dynamics of the model are described by a stochastic matrix — a nonnegative square matrix 𝑃 = 𝑃 [𝑖, 𝑗]
such that each row 𝑃 [𝑖, ⋅] sums to one.
The interpretation of 𝑃 is that 𝑃 [𝑖, 𝑗] represents the probability of transitioning from state 𝑖 to state 𝑗 in one unit of time.
In symbols,
ℙ{𝑋𝑡+1 = 𝑗 | 𝑋𝑡 = 𝑖} = 𝑃 [𝑖, 𝑗]
Equivalently,
• 𝑃 can be thought of as a family of distributions 𝑃 [𝑖, ⋅], one for each 𝑖 ∈ 𝑆
• 𝑃 [𝑖, ⋅] is the distribution of 𝑋𝑡+1 given 𝑋𝑡 = 𝑖
(As you probably recall, when using NumPy arrays, 𝑃 [𝑖, ⋅] is expressed as P[i,:])
In this section, we’ll allow 𝑆 to be a subset of ℝ, such as
• ℝ itself
• the positive reals (0, ∞)
• a bounded interval (𝑎, 𝑏)
The family of discrete distributions 𝑃 [𝑖, ⋅] will be replaced by a family of densities 𝑝(𝑥, ⋅), one for each 𝑥 ∈ 𝑆.
Analogous to the finite state case, 𝑝(𝑥, ⋅) is to be understood as the distribution (density) of 𝑋𝑡+1 given 𝑋𝑡 = 𝑥.
More formally, a stochastic kernel on 𝑆 is a function 𝑝 ∶ 𝑆 × 𝑆 → ℝ with the property that
1. 𝑝(𝑥, 𝑦) ≥ 0 for all 𝑥, 𝑦 ∈ 𝑆
2. ∫ 𝑝(𝑥, 𝑦)𝑑𝑦 = 1 for all 𝑥 ∈ 𝑆
1 (𝑦 − 𝑥)2
𝑝𝑤 (𝑥, 𝑦) ∶= √ exp {− } (2.1)
2𝜋 2
In the previous section, we made the connection between stochastic difference equation (2.2) and stochastic kernel (2.1).
In economics and time-series analysis we meet stochastic difference equations of all different shapes and sizes.
It will be useful for us if we have some systematic methods for converting stochastic difference equations into stochastic
kernels.
To this end, consider the generic (scalar) stochastic difference equation given by
This is a special case of (2.3) with 𝜇(𝑥) = 𝛼𝑥 and 𝜎(𝑥) = (𝛽 + 𝛾𝑥2 )1/2 .
Example 3: With stochastic production and a constant savings rate, the one-sector neoclassical growth model leads to a
law of motion for capital per worker such as
Here
• 𝑠 is the rate of savings
• 𝐴𝑡+1 is a production shock
1 𝑦 − 𝜇(𝑥)
𝑝(𝑥, 𝑦) = 𝜙( ) (2.7)
𝜎(𝑥) 𝜎(𝑥)
1 𝑦 − (1 − 𝛿)𝑥
𝑝(𝑥, 𝑦) = 𝜙( ) (2.8)
𝑠𝑓(𝑥) 𝑠𝑓(𝑥)
In this section of our lecture on finite Markov chains, we asked the following question: If
1. {𝑋𝑡 } is a Markov chain with stochastic matrix 𝑃
2. the distribution of 𝑋𝑡 is known to be 𝜓𝑡
then what is the distribution of 𝑋𝑡+1 ?
Letting 𝜓𝑡+1 denote the distribution of 𝑋𝑡+1 , the answer we gave was that
This intuitive equality states that the probability of being at 𝑗 tomorrow is the probability of visiting 𝑖 today and then
going on to 𝑗, summed over all possible 𝑖.
In the density case, we just replace the sum with an integral and probability mass functions with densities, yielding
Note: Unlike most operators, we write 𝑃 to the right of its argument, instead of to the left (i.e., 𝜓𝑃 instead of 𝑃 𝜓).
This is a common convention, with the intention being to maintain the parallel with the finite case — see here
With this notation, we can write (2.9) more succinctly as 𝜓𝑡+1 (𝑦) = (𝜓𝑡 𝑃 )(𝑦) for all 𝑦, or, dropping the 𝑦 and letting
“=” indicate equality of functions,
𝜓𝑡+1 = 𝜓𝑡 𝑃 (2.11)
Equation (2.11) tells us that if we specify a distribution for 𝜓0 , then the entire sequence of future distributions can be
obtained by iterating with 𝑃 .
It’s interesting to note that (2.11) is a deterministic difference equation.
Thus, by converting a stochastic difference equation such as (2.3) into a stochastic kernel 𝑝 and hence an operator 𝑃 , we
convert a stochastic difference equation into a deterministic one (albeit in a much higher dimensional space).
Note: Some people might be aware that discrete Markov chains are in fact a special case of the continuous Markov
chains we have just described. The reason is that probability mass functions are densities with respect to the counting
measure.
2.2.4 Computation
To learn about the dynamics of a given process, it’s useful to compute and study the sequences of densities generated by
the model.
One way to do this is to try to implement the iteration described by (2.10) and (2.11) using numerical integration.
However, to produce 𝜓𝑃 from 𝜓 via (2.10), you would need to integrate at every 𝑦, and there is a continuum of such 𝑦.
Another possibility is to discretize the model, but this introduces errors of unknown size.
A nicer alternative in the present setting is to combine simulation with an elegant estimator called the look-ahead estimator.
Let’s go over the ideas with reference to the growth model discussed above, the dynamics of which we repeat here for
convenience:
Our aim is to compute the sequence {𝜓𝑡 } associated with this model and fixed initial condition 𝜓0 .
To approximate 𝜓𝑡 by simulation, recall that, by definition, 𝜓𝑡 is the density of 𝑘𝑡 given 𝑘0 ∼ 𝜓0 .
If we wish to generate observations of this random variable, all we need to do is
1. draw 𝑘0 from the specified initial condition 𝜓0
1 𝑛
𝜓𝑡𝑛 (𝑦) = 𝑖
∑ 𝑝(𝑘𝑡−1 , 𝑦) (2.13)
𝑛 𝑖=1
1 𝑛 𝑖 𝑖
∑ 𝑝(𝑘𝑡−1 , 𝑦) → 𝔼𝑝(𝑘𝑡−1 , 𝑦) = ∫ 𝑝(𝑥, 𝑦)𝜓𝑡−1 (𝑥) 𝑑𝑥 = 𝜓𝑡 (𝑦)
𝑛 𝑖=1
2.2.5 Implementation
A class called LAE for estimating densities by this technique can be found in lae.py.
Given our use of the __call__ method, an instance of LAE acts as a callable object, which is essentially a function that
can store its own data (see this discussion).
This function returns the right-hand side of (2.13) using
• the data and stochastic kernel that it stores as its instance data
• the value 𝑦 as its argument
The function is vectorized, in the sense that if psi is such an instance and y is an array, then the call psi(y) acts
elementwise.
(This is the reason that we reshaped X and y inside the class — to make vectorization work)
Because the implementation is fully vectorized, it is about as efficient as it would be in C or Fortran.
2.2.6 Example
The following code is an example of usage for the stochastic growth model described above
# == Define parameters == #
s = 0.2
δ = 0.1
a_σ = 0.4 # A = exp(B) where B ~ N(0, a_σ)
α = 0.4 # We set f(k) = k**α
ψ_0 = beta(5, 5, scale=0.5) # Initial distribution
ϕ = lognorm(a_σ)
# == Generate T instances of LAE using this data, one for each date t == #
laes = [LAE(p, k[:, t]) for t in range(T)]
# == Plot == #
fig, ax = plt.subplots()
ygrid = np.linspace(0.01, 4.0, 200)
greys = [str(g) for g in np.linspace(0.0, 0.8, T)]
greys.reverse()
for ψ, g in zip(laes, greys):
ax.plot(ygrid, ψ(ygrid), color=g, lw=2, alpha=0.6)
ax.set_xlabel('capital')
ax.set_title(f'Density of $k_1$ (lighter) to $k_T$ (darker) for $T={T}$')
plt.show()
The figure shows part of the density sequence {𝜓𝑡 }, with each density computed via the look-ahead estimator.
Notice that the sequence of densities shown in the figure seems to be converging — more on this in just a moment.
Another quick comment is that each of these distributions could be interpreted as a cross-sectional distribution (recall
this discussion).
Up until now, we have focused exclusively on continuous state Markov chains where all conditional distributions 𝑝(𝑥, ⋅)
are densities.
As discussed above, not all distributions can be represented as densities.
If the conditional distribution of 𝑋𝑡+1 given 𝑋𝑡 = 𝑥 cannot be represented as a density for some 𝑥 ∈ 𝑆, then we need
a slightly different theory.
The ultimate option is to switch from densities to probability measures, but not all readers will be familiar with measure
theory.
We can, however, construct a fairly general theory using distribution functions.
To illustrate the issues, recall that Hopenhayn and Rogerson [Hopenhayn and Rogerson, 1993] study a model of firm
dynamics where individual firm productivity follows the exogenous process
IID
𝑋𝑡+1 = 𝑎 + 𝜌𝑋𝑡 + 𝜉𝑡+1 , where {𝜉𝑡 } ∼ 𝑁 (0, 𝜎2 )
If you think about it, you will see that for any given 𝑥 ∈ [0, 1], the conditional distribution of 𝑋𝑡+1 given 𝑋𝑡 = 𝑥 puts
positive probability mass on 0 and 1.
Hence it cannot be represented as a density.
What we can do instead is use cumulative distribution functions (cdfs).
To this end, set
This family of cdfs 𝐺(𝑥, ⋅) plays a role analogous to the stochastic kernel in the density case.
The distribution dynamics in (2.9) are then replaced by
Here 𝐹𝑡 and 𝐹𝑡+1 are cdfs representing the distribution of the current state and next period state.
The intuition behind (2.14) is essentially the same as for (2.9).
2.3.2 Computation
If you wish to compute these cdfs, you cannot use the look-ahead estimator as before.
Indeed, you should not use any density estimator, since the objects you are estimating/computing are not densities.
One good option is simulation as before, combined with the empirical distribution function.
2.4 Stability
In our lecture on finite Markov chains, we also studied stationarity, stability and ergodicity.
Here we will cover the same topics for the continuous case.
We will, however, treat only the density case (as in this section), where the stochastic kernel is a family of densities.
The general case is relatively similar — references are given below.
2.4. Stability 31
Advanced Quantitative Economics with Python
Analogous to the finite case, given a stochastic kernel 𝑝 and corresponding Markov operator as defined in (2.10), a density
𝜓∗ on 𝑆 is called stationary for 𝑃 if it is a fixed point of the operator 𝑃 .
In other words,
As with the finite case, if 𝜓∗ is stationary for 𝑃 , and the distribution of 𝑋0 is 𝜓∗ , then, in view of (2.11), 𝑋𝑡 will have
this same distribution for all 𝑡.
Hence 𝜓∗ is the stochastic equivalent of a steady state.
In the finite case, we learned that at least one stationary distribution exists, although there may be many.
When the state space is infinite, the situation is more complicated.
Even existence can fail very easily.
For example, the random walk model has no stationary density (see, e.g., EDTC, p. 210).
However, there are well-known conditions under which a stationary density 𝜓∗ exists.
With additional conditions, we can also get a unique stationary density (𝜓 ∈ 𝒟 and 𝜓 = 𝜓𝑃 ⟹ 𝜓 = 𝜓∗ ), and also
global convergence in the sense that
∀ 𝜓 ∈ 𝒟, 𝜓𝑃 𝑡 → 𝜓∗ as 𝑡 → ∞ (2.16)
This combination of existence, uniqueness and global convergence in the sense of (2.16) is often referred to as global
stability.
Under very similar conditions, we get ergodicity, which means that
1 𝑛
∑ ℎ(𝑋𝑡 ) → ∫ ℎ(𝑥)𝜓∗ (𝑥)𝑑𝑥 as 𝑛 → ∞ (2.17)
𝑛 𝑡=1
for any (measurable) function ℎ ∶ 𝑆 → ℝ such that the right-hand side is finite.
Note that the convergence in (2.17) does not depend on the distribution (or value) of 𝑋0 .
This is actually very important for simulation — it means we can learn about 𝜓∗ (i.e., approximate the right-hand side of
(2.17) via the left-hand side) without requiring any special knowledge about what to do with 𝑋0 .
So what are these conditions we require to get global stability and ergodicity?
In essence, it must be the case that
1. Probability mass does not drift off to the “edges” of the state space.
2. Sufficient “mixing” obtains.
For one such set of conditions see theorem 8.2.14 of EDTC.
In addition
• [Stokey et al., 1989] contains a classic (but slightly outdated) treatment of these topics.
• From the mathematical literature, [Lasota and MacKey, 1994] and [Meyn and Tweedie, 2009] give outstanding
in-depth treatments.
• Section 8.1.2 of EDTC provides detailed intuition, and section 8.3 gives additional references.
• EDTC, section 11.3.4 provides a specific treatment for the growth model we considered in this lecture.
As stated above, the growth model treated here is stable under mild conditions on the primitives.
• See EDTC, section 11.3.4 for more details.
We can see this stability in action — in particular, the convergence in (2.16) — by simulating the path of densities from
various initial conditions.
Here is such a figure.
All sequences are converging towards the same limit, regardless of their initial condition.
The details regarding initial conditions and so on are given in this exercise, where you are asked to replicate the figure.
In the preceding figure, each sequence of densities is converging towards the unique stationary density 𝜓∗ .
Even from this figure, we can get a fair idea what 𝜓∗ looks like, and where its mass is located.
However, there is a much more direct way to estimate the stationary density, and it involves only a slight modification of
the look-ahead estimator.
Let’s say that we have a model of the form (2.3) that is stable and ergodic.
Let 𝑝 be the corresponding stochastic kernel, as given in (2.7).
2.4. Stability 33
Advanced Quantitative Economics with Python
To approximate the stationary density 𝜓∗ , we can simply generate a long time-series 𝑋0 , 𝑋1 , … , 𝑋𝑛 and estimate 𝜓∗ via
1 𝑛
𝜓𝑛∗ (𝑦) = ∑ 𝑝(𝑋𝑡 , 𝑦) (2.18)
𝑛 𝑡=1
This is essentially the same as the look-ahead estimator (2.13), except that now the observations we generate are a single
time-series, rather than a cross-section.
The justification for (2.18) is that, with probability one as 𝑛 → ∞,
1 𝑛
∑ 𝑝(𝑋𝑡 , 𝑦) → ∫ 𝑝(𝑥, 𝑦)𝜓∗ (𝑥) 𝑑𝑥 = 𝜓∗ (𝑦)
𝑛 𝑡=1
where the convergence is by (2.17) and the equality on the right is by (2.15).
The right-hand side is exactly what we want to compute.
On top of this asymptotic result, it turns out that the rate of convergence for the look-ahead estimator is very good.
The first exercise helps illustrate this point.
2.5 Exercises
Exercise 2.5.1
Consider the simple threshold autoregressive model
IID
𝑋𝑡+1 = 𝜃|𝑋𝑡 | + (1 − 𝜃2 )1/2 𝜉𝑡+1 where {𝜉𝑡 } ∼ 𝑁 (0, 1) (2.19)
This is one of those rare nonlinear stochastic models where an analytical expression for the stationary density is available.
In particular, provided that |𝜃| < 1, there is a unique stationary density 𝜓∗ given by
𝜃𝑦
𝜓∗ (𝑦) = 2 𝜙(𝑦) Φ [ ] (2.20)
(1 − 𝜃2 )1/2
Here 𝜙 is the standard normal density and Φ is the standard normal cdf.
As an exercise, compute the look-ahead estimate of 𝜓∗ , as defined in (2.18), and compare it with 𝜓∗ in (2.20) to see
whether they are indeed close for large 𝑛.
In doing so, set 𝜃 = 0.8 and 𝑛 = 500.
The next figure shows the result of such a computation
The additional density (black line) is a nonparametric kernel density estimate, added to the solution for illustration.
(You can try to replicate it before looking at the solution if you want to)
As you can see, the look-ahead estimator is a much tighter fit than the kernel density estimator.
If you repeat the simulation you will see that this is consistently the case.
2.5. Exercises 35
Advanced Quantitative Economics with Python
ϕ = norm()
n = 500
θ = 0.8
# == Frequently used constants == #
d = np.sqrt(1 - θ**2)
δ = θ / d
def ψ_star(y):
"True stationary density of the TAR Model"
return 2 * norm.pdf(y) * norm.cdf(δ * y)
Z = ϕ.rvs(n)
X = np.empty(n)
for t in range(n-1):
X[t+1] = θ * np.abs(X[t]) + d * Z[t]
ψ_est = LAE(p, X)
k_est = gaussian_kde(X)
Exercise 2.5.2
Replicate the figure on global convergence shown above.
The densities come from the stochastic growth model treated at the start of the lecture.
Begin with the code found above.
Use the same parameters.
For the four initial distributions, use the shifted beta distributions
# == Define parameters == #
s = 0.2
δ = 0.1
a_σ = 0.4 # A = exp(B) where B ~ N(0, a_σ)
α = 0.4 # f(k) = k**α
ϕ = lognorm(a_σ)
(continues on next page)
2.5. Exercises 37
Advanced Quantitative Economics with Python
for i in range(4):
ax = axes[i]
ax.set_xlim(0, xmax)
ψ_0 = beta(5, 5, scale=0.5, loc=i*2) # Initial distribution
Exercise 2.5.3
A common way to compare distributions visually is with boxplots.
To illustrate, let’s generate three artificial data sets and compare them with a boxplot.
The three data sets we will use are:
{𝑋1 , … , 𝑋𝑛 } ∼ 𝐿𝑁 (0, 1), {𝑌1 , … , 𝑌𝑛 } ∼ 𝑁 (2, 1), and {𝑍1 , … , 𝑍𝑛 } ∼ 𝑁 (4, 1),
n = 500
x = np.random.randn(n) # N(0, 1)
x = np.exp(x) # Map x to lognormal
y = np.random.randn(n) + 2.0 # N(2, 1)
z = np.random.randn(n) + 4.0 # N(4, 1)
2.5. Exercises 39
Advanced Quantitative Economics with Python
Each data set is represented by a box, where the top and bottom of the box are the third and first quartiles of the data,
and the red line in the center is the median.
The boxes give some indication as to
• the location of probability mass for each sample
• whether the distribution is right-skewed (as is the lognormal distribution), etc
Now let’s put these ideas to use in a simulation.
Consider the threshold autoregressive model in (2.19).
We know that the distribution of 𝑋𝑡 will converge to (2.20) whenever |𝜃| < 1.
Let’s observe this convergence from different initial conditions using boxplots.
In particular, the exercise is to generate J boxplot figures, one for each initial condition 𝑋0 in
initial_conditions = np.linspace(8, 0, J)
Note the way we use vectorized code to simulate the 𝑘 time series for one boxplot all at once
n = 20
k = 5000
J = 8
θ = 0.9
d = np.sqrt(1 - θ**2)
δ = θ / d
for j in range(J):
axes[j].set_ylim(-4, 8)
axes[j].set_title(f'time series from t = {initial_conditions[j]}')
Z = np.random.randn(k, n)
X[:, 0] = initial_conditions[j]
for t in range(1, n):
X[:, t] = θ * np.abs(X[:, t-1]) + d * Z[:, t]
axes[j].boxplot(X)
plt.show()
2.5. Exercises 41
Advanced Quantitative Economics with Python
2.6 Appendix
2.6. Appendix 43
Advanced Quantitative Economics with Python
THREE
This lecture uses the Kalman filter to reformulate John F. Muth’s first paper [Muth, 1960] about rational expectations.
Muth used classical prediction methods to reverse engineer a stochastic process that renders optimal Milton Friedman’s
[Friedman, 1956] “adaptive expectations” scheme.
Milton Friedman [Friedman, 1956] (1956) posited that consumer’s forecast their future disposable income with the adap-
tive expectations scheme
∞
∗
𝑦𝑡+𝑖,𝑡 = 𝐾 ∑(1 − 𝐾)𝑗 𝑦𝑡−𝑗 (3.1)
𝑗=0
∗
where 𝐾 ∈ (0, 1) and 𝑦𝑡+𝑖,𝑡 is a forecast of future 𝑦 over horizon 𝑖.
Milton Friedman justified the exponential smoothing forecasting scheme (3.1) informally, noting that it seemed a plau-
sible way to use past income to forecast future income.
In his first paper about rational expectations, John F. Muth [Muth, 1960] reverse-engineered a univariate stochastic
∞
process {𝑦𝑡 }𝑡=−∞ for which Milton Friedman’s adaptive expectations scheme gives linear least forecasts of 𝑦𝑡+𝑗 for any
horizon 𝑖.
Muth sought a setting and a sense in which Friedman’s forecasting scheme is optimal.
That is, Muth asked for what optimal forecasting question is Milton Friedman’s adaptive expectation scheme the answer.
Muth (1960) used classical prediction methods based on lag-operators and 𝑧-transforms to find the answer to his question.
Please see lectures Classical Control with Linear Algebra and Classical Filtering and Prediction with Linear Algebra for an
introduction to the classical tools that Muth used.
45
Advanced Quantitative Economics with Python
Rather than using those classical tools, in this lecture we apply the Kalman filter to express the heart of Muth’s analysis
concisely.
The lecture First Look at Kalman Filter describes the Kalman filter.
We’ll use limiting versions of the Kalman filter corresponding to what are called stationary values in that lecture.
Suppose that an observable 𝑦𝑡 is the sum of an unobserved random walk 𝑥𝑡 and an IID shock 𝜖2,𝑡 :
𝑥𝑡+1 = 𝑥𝑡 + 𝜎𝑥 𝜖1,𝑡+1
(3.2)
𝑦𝑡 = 𝑥𝑡 + 𝜎𝑦 𝜖2,𝑡
where
𝜖
[ 1,𝑡+1 ] ∼ 𝒩(0, 𝐼)
𝜖2,𝑡
is an IID process.
Note: A property of the state-space representation (3.2) is that in general neither 𝜖1,𝑡 nor 𝜖2,𝑡 is in the space spanned by
square-summable linear combinations of 𝑦𝑡 , 𝑦𝑡−1 , ….
𝜖
In general [ 1,𝑡 ] has more information about future 𝑦𝑡+𝑗 ’s than is contained in 𝑦𝑡 , 𝑦𝑡−1 , ….
𝜖2𝑡
We can use the asymptotic or stationary values of the Kalman gain and the one-step-ahead conditional state covariance
matrix to compute a time-invariant innovations representation
𝑥𝑡+1
̂ = 𝑥𝑡̂ + 𝐾𝑎𝑡
(3.3)
𝑦𝑡 = 𝑥𝑡̂ + 𝑎𝑡
where 𝑥𝑡̂ = 𝐸[𝑥𝑡 |𝑦𝑡−1 , 𝑦𝑡−2 , …] and 𝑎𝑡 = 𝑦𝑡 − 𝐸[𝑦𝑡 |𝑦𝑡−1 , 𝑦𝑡−2 , …].
Note: A key property about an innovations representation is that 𝑎𝑡 is in the space spanned by square summable linear
combinations of 𝑦𝑡 , 𝑦𝑡−1 , ….
For more ramifications of this property, see the lectures Shock Non-Invertibility and Recursive Models of Dynamic Linear
Economies.
Later we’ll stack these state-space systems (3.2) and (3.3) to display some classic findings of Muth.
But first, let’s create an instance of the state-space system (3.2) then apply the quantecon Kalman class, then uses it to
construct the associated “innovations representation”
Now we want to map the time-invariant innovations representation (3.3) and the original state-space system (3.2) into a
convenient form for deducing the impulse responses from the original shocks to the 𝑥𝑡 and 𝑥𝑡̂ .
Putting both of these representations into a single state-space system is yet another application of the insight that “finding
the state is an art”.
We’ll define a state vector and appropriate state-space matrices that allow us to represent both systems in one fell swoop.
Note that
𝑎𝑡 = 𝑥𝑡 + 𝜎𝑦 𝜖2,𝑡 − 𝑥𝑡̂
so that
𝑥𝑡+1
̂ = 𝑥𝑡̂ + 𝐾(𝑥𝑡 + 𝜎𝑦 𝜖2,𝑡 − 𝑥𝑡̂ )
= (1 − 𝐾)𝑥𝑡̂ + 𝐾𝑥𝑡 + 𝐾𝜎𝑦 𝜖2,𝑡
𝑥𝑡+1 1 0 0 𝑥𝑡 𝜎𝑥 0
⎤ = ⎡𝐾 𝜖1,𝑡+1
⎡ 𝑥̂
⎢ 𝑡+1 ⎥ ⎢ (1 − 𝐾) 𝐾𝜎𝑦 ⎤ ⎡ 𝑥̂ ⎤ + ⎡ 0
⎥⎢ 𝑡 ⎥ ⎢ 0⎤
⎥ [𝜖 ]
⎣𝜖2,𝑡+1 ⎦ ⎣ 0 0 0 ⎦ ⎣𝜖2,𝑡 ⎦ ⎣ 0 1⎦ 2,𝑡+1
𝑥
𝑦 1 0 𝜎𝑦 ⎡ 𝑡 ⎤
[ 𝑡] = [ ] ⎢ 𝑥𝑡̂ ⎥
𝑎𝑡 1 −1 𝜎𝑦
⎣𝜖2,𝑡 ⎦
𝜖
is a state-space system that tells us how the shocks [ 1,𝑡+1 ] affect states 𝑥𝑡+1
̂ , 𝑥𝑡 , the observable 𝑦𝑡 , and the innovation
𝜖2,𝑡+1
𝑎𝑡 .
With this tool at our disposal, let’s form the composite system and simulate it
Now that we have simulated our joint system, we have 𝑥𝑡 , 𝑥𝑡̂ , and 𝑦𝑡 .
We can now investigate how these variables are related by plotting some key objects.
First, let’s plot the hidden state 𝑥𝑡 and the filtered version 𝑥𝑡̂ that is linear-least squares projection of 𝑥𝑡 on the history
𝑦𝑡−1 , 𝑦𝑡−2 , …
fig, ax = plt.subplots()
ax.plot(xf[0, :], label="$x_t$")
ax.plot(xf[1, :], label="Filtered $x_t$")
ax.legend()
ax.set_xlabel("Time")
ax.set_title(r"$x$ vs $\hat{x}$")
plt.show()
fig, ax = plt.subplots()
ax.plot(yf[0, :], label="y")
ax.plot(xf[0, :], label="x")
ax.legend()
ax.set_title(r"$x$ and $y$")
ax.set_xlabel("Time")
plt.show()
We see above that 𝑦 seems to look like white noise around the values of 𝑥.
3.5.1 Innovations
Recall that we wrote down the innovation representation that depended on 𝑎𝑡 . We now plot the innovations {𝑎𝑡 }:
fig, ax = plt.subplots()
ax.plot(yf[1, :], label="a")
ax.legend()
ax.set_title(r"Innovation $a_t$")
ax.set_xlabel("Time")
plt.show()
fig, ax = plt.subplots(2)
ax[0].plot(coefs_ma_array, label="MA")
ax[0].legend()
ax[1].plot(coefs_var_array, label="VAR")
(continues on next page)
plt.show()
The moving average coefficients in the top panel show tell-tale signs of 𝑦𝑡 being a process whose first difference is a
first-order autoregression.
The autoregressive coefficients decline geometrically with decay rate (1 − 𝐾).
These are exactly the target outcomes that Muth (1960) aimed to reverse engineer
FOUR
In addition to what’s in Anaconda, this lecture will need the following libraries:
4.1 Overview
In this lecture we discuss a family of dynamic programming problems with the following features:
1. a discrete state space and discrete choices (actions)
2. an infinite horizon
3. discounted rewards
4. Markov state transitions
We call such problems discrete dynamic programs or discrete DPs.
Discrete DPs are the workhorses in much of modern quantitative economics, including
• monetary economics
• search and labor economics
• household savings and consumption theory
• investment theory
• asset pricing
• industrial organization, etc.
When a given model is not inherently discrete, it is common to replace it with a discretized version in order to use discrete
DP techniques.
This lecture covers
• the theory of dynamic programming in a discrete setting, plus examples and applications
• a powerful set of routines for solving discrete DPs from the QuantEcon code library
Let’s start with some imports:
import numpy as np
import matplotlib.pyplot as plt
import quantecon as qe
import scipy.sparse as sparse
(continues on next page)
53
Advanced Quantitative Economics with Python
4.1.2 Code
4.1.3 References
For background reading on dynamic programming and additional applications, see, for example,
• [Ljungqvist and Sargent, 2018]
• [Hernandez-Lerma and Lasserre, 1996], section 3.5
• [Puterman, 2005]
• [Stokey et al., 1989]
• [Rust, 1996]
• [Miranda and Fackler, 2002]
• EDTC, chapter 5
Loosely speaking, a discrete DP is a maximization problem with an objective function of the form
∞
𝔼 ∑ 𝛽 𝑡 𝑟(𝑠𝑡 , 𝑎𝑡 ) (4.1)
𝑡=0
where
• 𝑠𝑡 is the state variable
• 𝑎𝑡 is the action
• 𝛽 is a discount factor
• 𝑟(𝑠𝑡 , 𝑎𝑡 ) is interpreted as a current reward when the state is 𝑠𝑡 and the action chosen is 𝑎𝑡
Each pair (𝑠𝑡 , 𝑎𝑡 ) pins down transition probabilities 𝑄(𝑠𝑡 , 𝑎𝑡 , 𝑠𝑡+1 ) for the next period state 𝑠𝑡+1 .
Thus, actions influence not only current rewards but also the future time path of the state.
The essence of dynamic programming problems is to trade off current rewards vs favorable positioning of the future state
(modulo randomness).
Examples:
• consuming today vs saving and accumulating assets
• accepting a job offer today vs seeking a better one in the future
• exercising an option now vs waiting
4.2.1 Policies
The most fruitful way to think about solutions to discrete DP problems is to compare policies.
In general, a policy is a randomized map from past actions and states to current action.
In the setting formalized below, it suffices to consider so-called stationary Markov policies, which consider only the current
state.
In particular, a stationary Markov policy is a map 𝜎 from states to actions
• 𝑎𝑡 = 𝜎(𝑠𝑡 ) indicates that 𝑎𝑡 is the action to be taken in state 𝑠𝑡
It is known that, for any arbitrary policy, there exists a stationary Markov policy that dominates it at least weakly.
• See section 5.5 of [Puterman, 2005] for discussion and proofs.
In what follows, stationary Markov policies are referred to simply as policies.
The aim is to find an optimal policy, in the sense of one that maximizes (4.1).
Let’s now step through these ideas more carefully.
SA ∶= {(𝑠, 𝑎) ∣ 𝑠 ∈ 𝑆, 𝑎 ∈ 𝐴(𝑠)}
3. A reward function 𝑟 ∶ SA → ℝ.
4. A transition probability function 𝑄 ∶ SA → Δ(𝑆), where Δ(𝑆) is the set of probability distributions over 𝑆.
5. A discount factor 𝛽 ∈ [0, 1).
We also use the notation 𝐴 ∶= ⋃𝑠∈𝑆 𝐴(𝑠) = {0, … , 𝑚 − 1} and call this set the action space.
A policy is a function 𝜎 ∶ 𝑆 → 𝐴.
A policy is called feasible if it satisfies 𝜎(𝑠) ∈ 𝐴(𝑠) for all 𝑠 ∈ 𝑆.
Denote the set of all feasible policies by Σ.
If a decision-maker uses a policy 𝜎 ∈ Σ, then
• the current reward at time 𝑡 is 𝑟(𝑠𝑡 , 𝜎(𝑠𝑡 ))
• the probability that 𝑠𝑡+1 = 𝑠′ is 𝑄(𝑠𝑡 , 𝜎(𝑠𝑡 ), 𝑠′ )
For each 𝜎 ∈ Σ, define
• 𝑟𝜎 by 𝑟𝜎 (𝑠) ∶= 𝑟(𝑠, 𝜎(𝑠)))
• 𝑄𝜎 by 𝑄𝜎 (𝑠, 𝑠′ ) ∶= 𝑄(𝑠, 𝜎(𝑠), 𝑠′ )
Notice that 𝑄𝜎 is a stochastic matrix on 𝑆.
It gives transition probabilities of the controlled chain when we follow policy 𝜎.
If we think of 𝑟𝜎 as a column vector, then so is 𝑄𝑡𝜎 𝑟𝜎 , and the 𝑠-th row of the latter has the interpretation
Comments
• {𝑠𝑡 } ∼ 𝑄𝜎 means that the state is generated by stochastic matrix 𝑄𝜎 .
• See this discussion on computing expectations of Markov chains for an explanation of the expression in (4.2).
Notice that we’re not really distinguishing between functions from 𝑆 to ℝ and vectors in ℝ𝑛 .
This is natural because they are in one to one correspondence.
Let 𝑣𝜎 (𝑠) denote the discounted sum of expected reward flows from policy 𝜎 when the initial state is 𝑠.
To calculate this quantity we pass the expectation through the sum in (4.1) and use (4.2) to get
∞
𝑣𝜎 (𝑠) = ∑ 𝛽 𝑡 (𝑄𝑡𝜎 𝑟𝜎 )(𝑠) (𝑠 ∈ 𝑆)
𝑡=0
This function is called the policy value function for the policy 𝜎.
The optimal value function, or simply value function, is the function 𝑣∗ ∶ 𝑆 → ℝ defined by
(We can use max rather than sup here because the domain is a finite set)
A policy 𝜎 ∈ Σ is called optimal if 𝑣𝜎 (𝑠) = 𝑣∗ (𝑠) for all 𝑠 ∈ 𝑆.
Given any 𝑤 ∶ 𝑆 → ℝ, a policy 𝜎 ∈ Σ is called 𝑤-greedy if
As discussed in detail below, optimal policies are precisely those that are 𝑣∗ -greedy.
Now that the theory has been set out, let’s turn to solution methods.
The code for solving discrete DPs is available in ddp.py from the QuantEcon.py code library.
It implements the three most important solution methods for discrete dynamic programs, namely
• value function iteration
• policy function iteration
• modified policy function iteration
Let’s briefly review these algorithms and their implementation.
Perhaps the most familiar method for solving all manner of dynamic programs is value function iteration.
This algorithm uses the fact that the Bellman operator 𝑇 is a contraction mapping with fixed point 𝑣∗ .
Hence, iterative application of 𝑇 to any initial function 𝑣0 ∶ 𝑆 → ℝ converges to 𝑣∗ .
The details of the algorithm can be found in the appendix.
This routine, also known as Howard’s policy improvement algorithm, exploits more closely the particular structure of a
discrete DP problem.
Each iteration consists of
1. A policy evaluation step that computes the value 𝑣𝜎 of a policy 𝜎 by solving the linear equation 𝑣 = 𝑇𝜎 𝑣.
2. A policy improvement step that computes a 𝑣𝜎 -greedy policy.
In the current setting, policy iteration computes an exact optimal policy in finitely many iterations.
• See theorem 10.2.6 of EDTC for a proof.
The details of the algorithm can be found in the appendix.
Modified policy iteration replaces the policy evaluation step in policy iteration with “partial policy evaluation”.
The latter computes an approximation to the value of a policy 𝜎 by iterating 𝑇𝜎 for a specified number of times.
This approach can be useful when the state space is very large and the linear system in the policy evaluation step of policy
iteration is correspondingly difficult to solve.
The details of the algorithm can be found in the appendix.
𝑠′ = 𝑎 + 𝑈 where 𝑈 ∼ 𝑈 [0, … , 𝐵]
This information will be used to create an instance of DiscreteDP by passing the following information
1. An 𝑛 × 𝑚 reward array 𝑅.
2. An 𝑛 × 𝑚 × 𝑛 transition probability array 𝑄.
3. A discount factor 𝛽.
For 𝑅 we set 𝑅[𝑠, 𝑎] = 𝑢(𝑠 − 𝑎) if 𝑎 ≤ 𝑠 and −∞ otherwise.
For 𝑄 we follow the rule in (4.3).
Note:
• The feasibility constraint is embedded into 𝑅 by setting 𝑅[𝑠, 𝑎] = −∞ for 𝑎 ∉ 𝐴(𝑠).
• Probability distributions for (𝑠, 𝑎) with 𝑎 ∉ 𝐴(𝑠) can be arbitrary.
class SimpleOG:
self.populate_Q()
self.populate_R()
def populate_R(self):
"""
Populate the R matrix, with R[s, a] = -np.inf for infeasible
state-action pairs.
"""
for s in range(self.n):
for a in range(self.m):
self.R[s, a] = self.u(s - a) if a <= s else -np.inf
def populate_Q(self):
"""
Populate the Q matrix by setting
for a in range(self.m):
self.Q[:, a, a:(a + self.B + 1)] = 1.0 / (self.B + 1)
results = ddp.solve(method='policy_iteration')
dir(results)
(In IPython version 4.0 and above you can also type results. and hit the tab key)
The most important attributes are v, the value function, and σ, the optimal policy
results.v
results.sigma
array([0, 0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 4, 5, 5, 5, 5])
Since we’ve used policy iteration, these results will be exact unless we hit the iteration bound max_iter.
Let’s make sure this didn’t happen
results.max_iter
250
results.num_iter
Another interesting object is results.mc, which is the controlled chain defined by 𝑄𝜎∗ , where 𝜎∗ is the optimal policy.
In other words, it gives the dynamics of the state when the agent follows the optimal policy.
Since this object is an instance of MarkovChain from QuantEcon.py (see this lecture for more discussion), we can easily
simulate it, compute its stationary distribution and so on.
results.mc.stationary_distributions
If we look at the bar graph we can see the rightward shift in probability mass
def u(c):
return c**α
s_indices = []
a_indices = []
Q = []
R = []
b = 1.0 / (B + 1)
for s in range(n):
for a in range(min(M, s) + 1): # All feasible a at this s
s_indices.append(s)
a_indices.append(a)
q = np.zeros(n)
q[a:(a + B + 1)] = b # b on these values, otherwise 0
Q.append(q)
R.append(u(s - a))
For larger problems, you might need to write this code more efficiently by vectorizing or using Numba.
4.5 Exercises
In the stochastic optimal growth lecture from our introductory lecture series, we solve a benchmark model that has an
analytical solution.
The exercise is to replicate this solution using DiscreteDP.
4.5. Exercises 63
Advanced Quantitative Economics with Python
4.6 Solutions
4.6.1 Setup
α = 0.65
f = lambda k: k**α
u = np.log
β = 0.95
Here we want to solve a finite state version of the continuous state model above.
We discretize the state space into a grid of size grid_size=500, from 10−6 to grid_max=2
grid_max = 2
grid_size = 500
grid = np.linspace(1e-6, grid_max, grid_size)
We choose the action to be the amount of capital to save for the next period (the state is the capital stock at the beginning
of the period).
Thus the state indices and the action indices are both 0, …, grid_size-1.
Action (indexed by) a is feasible at state (indexed by) s if and only if grid[a] < f([grid[s]) (zero consumption
is not allowed because of the log utility).
Thus the Bellman equation is:
# State-action indices
s_indices, a_indices = np.where(C > 0)
print(L)
print(s_indices)
print(a_indices)
118841
[ 0 1 1 ... 499 499 499]
[ 0 0 1 ... 389 390 391]
R = u(C[s_indices, a_indices])
(Degenerate) transition probability matrix Q (of shape (L, grid_size)), where we choose the scipy.sparse.lil_matrix
format, while any format will do (internally it will be converted to the csr format):
Q = sparse.lil_matrix((L, grid_size))
Q[np.arange(L), a_indices] = 1
(If you are familiar with the data structure of scipy.sparse.csr_matrix, the following is the most efficient way to create the
Q matrix in the current case)
# data = np.ones(L)
# indptr = np.arange(L+1)
# Q = sparse.csr_matrix((data, a_indices, indptr), shape=(L, grid_size))
Notes
Here we intensively vectorized the operations on arrays to simplify the code.
As noted, however, vectorization is memory consumptive, and it can be prohibitively so for grids with large size.
res = ddp.solve(method='policy_iteration')
v, σ, num_iter = res.v, res.sigma, res.num_iter
num_iter
10
Note that sigma contains the indices of the optimal capital stocks to save for the next period. The following translates
sigma to the corresponding consumption vector.
def v_star(k):
return c1 + c2 * np.log(k)
def c_star(k):
return (1 - ab) * k**α
Let us compare the solution of the discrete model with that of the original continuous model
4.6. Solutions 65
Advanced Quantitative Economics with Python
np.abs(v - v_star(grid)).max()
121.49819147053378
np.abs(v - v_star(grid))[1:].max()
0.012681735127500815
np.abs(c - c_star(grid)).max()
0.003826523100010082
In fact, the optimal consumption obtained in the discrete version is not really monotone, but the decrements are quite
small:
diff = np.diff(c)
(diff >= 0).all()
False
174
np.abs(diff[dec_ind]).max()
0.001961853339766839
True
Value Iteration
ddp.epsilon = 1e-4
ddp.max_iter = 500
res1 = ddp.solve(method='value_iteration')
res1.num_iter
294
np.array_equal(σ, res1.sigma)
True
4.6. Solutions 67
Advanced Quantitative Economics with Python
res2 = ddp.solve(method='modified_policy_iteration')
res2.num_iter
16
np.array_equal(σ, res2.sigma)
True
Speed Comparison
%timeit ddp.solve(method='value_iteration')
%timeit ddp.solve(method='policy_iteration')
%timeit ddp.solve(method='modified_policy_iteration')
94.9 ms ± 360 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
9.34 ms ± 16.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
11.3 ms ± 59.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
As is often the case, policy iteration and modified policy iteration are much faster than value iteration.
Let us first visualize the convergence of the value iteration algorithm as in the lecture, where we use ddp.
bellman_operator implemented as a method of DiscreteDP
plt.show()
We next plot the consumption policies along with the value iteration
4.6. Solutions 69
Advanced Quantitative Economics with Python
/home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages/quantecon/_
↪compute_fp.py:152: RuntimeWarning: max_iter attained before convergence in␣
↪compute_fixed_point
warnings.warn(_non_convergence_msg, RuntimeWarning)
4.6. Solutions 71
Advanced Quantitative Economics with Python
Finally, let us work on Exercise 2, where we plot the trajectories of the capital stock for three different discount factors,
0.9, 0.94, and 0.98, with initial condition 𝑘0 = 0.1.
sample_size = 25
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel("time")
ax.set_ylabel("capital")
ax.set_ylim(0.10, 0.30)
ax.legend(loc='lower right')
plt.show()
This appendix covers the details of the solution algorithms implemented for DiscreteDP.
We will make use of the following notions of approximate optimality:
• For 𝜀 > 0, 𝑣 is called an 𝜀-approximation of 𝑣∗ if ‖𝑣 − 𝑣∗ ‖ < 𝜀.
• A policy 𝜎 ∈ Σ is called 𝜀-optimal if 𝑣𝜎 is an 𝜀-approximation of 𝑣∗ .
The DiscreteDP value iteration method implements value function iteration as follows
1. Choose any 𝑣0 ∈ ℝ𝑛 , and specify 𝜀 > 0; set 𝑖 = 0.
2. Compute 𝑣𝑖+1 = 𝑇 𝑣𝑖 .
3. If ‖𝑣𝑖+1 − 𝑣𝑖 ‖ < [(1 − 𝛽)/(2𝛽)]𝜀, then go to step 4; otherwise, set 𝑖 = 𝑖 + 1 and go to step 2.
4. Compute a 𝑣𝑖+1 -greedy policy 𝜎, and return 𝑣𝑖+1 and 𝜎.
Given 𝜀 > 0, the value iteration algorithm
• terminates in a finite number of iterations
• returns an 𝜀/2-approximation of the optimal value function and an 𝜀-optimal policy function (unless iter_max
is reached)
(While not explicit, in the actual implementation each algorithm is terminated if the number of iterations reaches
iter_max)
LQ Control
75
CHAPTER
FIVE
5.1 Overview
In the linear-quadratic permanent income of consumption smoothing model described in this quantecon lecture, a scalar
parameter 𝛽 ∈ (0, 1) plays two roles:
• it is a discount factor that the consumer applies to future utilities from consumption
• it is the reciprocal of the gross interest rate on risk-free one-period loans
That 𝛽 plays these two roles is essential in delivering the outcome that, regardless of the stochastic process that describes
his non-financial income, the consumer chooses to make consumption follow a random walk (see [Hall, 1978]).
In this lecture, we assign a third role to 𝛽:
• it describes a first-order moving average process for the growth in non-financial income
We study two consumers who have exactly the same nonfinancial income process and who both conform to the linear-
quadratic permanent income of consumption smoothing model described here.
The two consumers have different information about their future nonfinancial incomes.
A better informed consumer each period receives news in the form of a shock that simultaneously affects both today’s
nonfinancial income and the present value of future nonfinancial incomes in a particular way.
A less informed consumer each period receives a shock that equals the part of today’s nonfinancial income that could not
be forecast from past values of nonfinancial income.
Even though they receive exactly the same nonfinancial incomes each period, our two consumers behave differently
because they have different information about their future nonfinancial incomes.
The second consumer receives less information about future nonfinancial incomes in a sense that we shall make precise.
This difference in their information sets manifests itself in their responding differently to what they regard as time 𝑡
information shocks.
Thus, although at each date they receive exactly the same histories of nonfinancial income, our two consumers receive
different shocks or news about their future nonfinancial incomes.
77
Advanced Quantitative Economics with Python
We study consequences of endowing a consumer with one of two alternative representations for the change in the con-
sumer’s nonfinancial income 𝑦𝑡+1 − 𝑦𝑡 .
For both types of consumer, a parameter 𝛽 ∈ (0, 1) plays three roles.
It appears
• as a discount factor applied to future expected one-period utilities,
• as the reciprocal of a gross interest rate on one-period loans, and
• as a parameter in a first-order moving average that equals the increment in a consumer’s non-financial income
The first representation, which we shall sometimes refer to as the more informative representation, is
where {𝜖𝑡 } is an i.i.d. normally distributed scalar process with means of zero and contemporaneous variances 𝜎𝜖2 .
This representation of the process is used by a consumer who at time 𝑡 knows both 𝑦𝑡 and the shock 𝜖𝑡 and can use both
of them to forecast future 𝑦𝑡+𝑗 ’s.
As we’ll see below, representation (5.1) has the peculiar property that a positive shock 𝜖𝑡+1 leaves the discounted present
value of the consumer’s financial income at time 𝑡 + 1 unaltered.
The second representation of the same {𝑦𝑡 } process is
where {𝑎𝑡 } is another i.i.d. normally distributed scalar process, with means of zero and now variances 𝜎𝑎2 > 𝜎𝜖2 .
The i.i.d. shock variances are related by
so that the variance of the innovation exceeds the variance of the original shock by a multiplicative factor 𝛽 −2 .
Representation (5.2) is the innovations representation of equation (5.1) associated with Kalman filtering theory.
To see how this works, note that equating representations (5.1) and (5.2) for 𝑦𝑡+1 −𝑦𝑡 implies 𝜖𝑡+1 −𝛽 −1 𝜖𝑡 = 𝑎𝑡+1 −𝛽𝑎𝑡 ,
which in turn implies
Solving this difference equation backwards for 𝑎𝑡+1 gives, after a few lines of algebra,
∞
𝑎𝑡+1 = 𝜖𝑡+1 + (𝛽 − 𝛽 −1 ) ∑ 𝛽 𝑗 𝜖𝑡−𝑗 (5.3)
𝑗=0
∞
where 𝐿 is the one-period lag operator, ℎ(𝐿) = ∑𝑗=0 ℎ𝑗 𝐿𝑗 , 𝐼 is the identity operator, and
𝐼 − 𝛽 −1 𝐿
ℎ(𝐿) =
𝐼 − 𝛽𝐿
Let 𝑔𝑗 ≡ 𝐸𝑧𝑡 𝑧𝑡−𝑗 be the 𝑗th autocovariance of the {𝑦𝑡 − 𝑦𝑡−1 } process.
Using calculations in the quantecon lecture, where 𝑧 ∈ 𝐶 is a complex variable, the covariance generating function
∞
𝑔(𝑧) = ∑𝑗=−∞ 𝑔𝑗 𝑧 𝑗 of the {𝑦𝑡 − 𝑦𝑡−1 } process equals
𝜎𝑎2 = 𝛽 −1 𝜎𝜖2 .
To verify these claims, just notice that 𝑔(𝑧) = 𝛽 −2 𝜎𝜖2 implies that
• 𝑔0 = 𝛽 −2 𝜎𝜖2 , and
• 𝑔𝑗 = 0 for 𝑗 ≠ 0.
Alternatively, if you are uncomfortable with covariance generating functions, note that we can directly calculate 𝜎𝑎2 from
formula (5.3) according to
∞
𝜎𝑎2 = 𝜎𝜖2 + [1 + (𝛽 − 𝛽 −1 )2 ∑ 𝛽 2𝑗 ] = 𝛽 −1 𝜎𝜖2 .
𝑗=0
We can also use the the Kalman filter to obtain representation (5.2) from representation (5.1).
Thus, from equations associated with the Kalman filter, it can be verified that the steady-state Kalman gain 𝐾 = 𝛽 2 and
the steady state conditional covariance
In a little more detail, let 𝑧𝑡 = 𝑦𝑡 − 𝑦𝑡−1 and form the state-space representation
𝜖𝑡+1
̂ = 0𝜖𝑡̂ + 𝐾𝑎𝑡+1
𝑧𝑡+1 = −𝛽𝑎𝑡 + 𝑎𝑡+1
By applying formulas for the steady-state Kalman filter, by hand it is possible to verify that 𝐾 = 𝛽 2 , 𝜎𝑎2 = 𝛽 −2 𝜎𝜖2 = 𝛽 −2 ,
and Σ = (1 − 𝛽 2 )𝜎𝜖2 .
Alternatively, we can obtain these formulas via the classical filtering theory described in this lecture.
Representation (5.1) is cast in terms of a news shock 𝜖𝑡+1 that represents a shock to nonfinancial income coming from
taxes, transfers, and other random sources of income changes known to a well-informed person who perhaps has all sorts
of information about the income process.
Representation (5.2) for the same income process is driven by shocks 𝑎𝑡 that contain less information than the news shock
𝜖𝑡 .
Representation (5.2) is called the innovations representation for the {𝑦𝑡 − 𝑦𝑡−1 } process.
It is cast in terms of what time series statisticians call the innovation or fundamental shock that emerges from apply-
ing the theory of optimally predicting nonfinancial income based solely on the information in past levels of growth in
nonfinancial income.
Fundamental for the 𝑦𝑡 process means that the shock 𝑎𝑡 can be expressed as a square-summable linear combination of
𝑦𝑡 , 𝑦𝑡−1 , ….
The shock 𝜖𝑡 is not fundamental because it has more information about the future of the {𝑦𝑡 − 𝑦𝑡−1 } process than is
contained in 𝑎𝑡 .
Representation (5.3) reveals the important fact that the original shock 𝜖𝑡 contains more information about future 𝑦’s than
is contained in the semi-infinite history 𝑦𝑡 = [𝑦𝑡 , 𝑦𝑡−1 , …].
Staring at representation (5.3) for 𝑎𝑡+1 shows that it consists both of new news 𝜖𝑡+1 as well as a long moving average
∞
(𝛽 − 𝛽 −1 ) ∑𝑗=0 𝛽 𝑗 𝜖𝑡−𝑗 of old news.
The more information representation (5.1) asserts that a shock 𝜖𝑡 results in an impulse response to nonfinancial income
of 𝜖𝑡 times the sequence
1, 1 − 𝛽 −1 , 1 − 𝛽 −1 , …
so that a shock that increases nonfinancial income 𝑦𝑡 by 𝜖𝑡 at time 𝑡 is followed by a change in future 𝑦 of 𝜖𝑡 times
1 − 𝛽 −1 < 0 in all subsequent periods.
Because 1 − 𝛽 −1 < 0, this means that a positive shock of 𝜖𝑡 today raises income at time 𝑡 by 𝜖𝑡 and then permanently
decreases all future incomes by (𝛽 −1 − 1)𝜖𝑡 .
This pattern precisely describes the following mental experiment:
• The consumer receives a government transfer of 𝜖𝑡 at time 𝑡.
• The government finances the transfer by issuing a one-period bond on which it pays a gross one-period risk-free
interest rate equal to 𝛽 −1 .
• In each future period, the government rolls over the one-period bond and so continues to borrow 𝜖𝑡 forever.
• The government imposes a lump-sum tax on the consumer in order to pay just the current interest on the original
bond and its rolled over successors.
• Thus, in periods 𝑡 + 1, 𝑡 + 2, …, the government levies a lump-sum tax on the consumer of 𝛽 −1 − 1 that is just
enough to pay the interest on the bond.
0
The present value of the impulse response or moving average coefficients equals 𝑑𝜖 (𝐿) = 1−𝛽 = 0, a fact that we’ll see
again below.
Representation (5.2), i.e., the innovations representation, asserts that a shock 𝑎𝑡 results in an impulse response to nonfi-
nancial income of 𝑎𝑡 times
1, 1 − 𝛽, 1 − 𝛽, …
so that a shock that increases income 𝑦𝑡 by 𝑎𝑡 at time 𝑡 can be expected to be followed by an increase in 𝑦𝑡+𝑗 of 𝑎𝑡 times
1 − 𝛽 > 0 in all future periods 𝑗 = 1, 2, ….
1−𝛽2
The present value of the impulse response or moving average coefficients for representation (5.2) is 𝑑𝑎 (𝛽) = 1−𝛽 =
(1 + 𝛽), another fact that will be important below.
Notice that reprentation (5.1), namely, 𝑦𝑡+1 − 𝑦𝑡 = −𝛽 −1 𝜖𝑡 + 𝜖𝑡+1 implies the linear difference equation
𝜖𝑡 = 𝛽𝜖𝑡+1 − 𝛽(𝑦𝑡+1 − 𝑦𝑡 ).
This equation shows that 𝜖𝑡 equals 𝛽 times the one-step-backwards error in optimally backcasting 𝑦𝑡 based on the semi-
𝑡
infinite future 𝑦+ ≡ [𝑦𝑡+1 , 𝑦𝑡+2 , …] via the optimal backcasting formula
∞
𝑡
𝐸[𝑦𝑡 |𝑦+ ] = (1 − 𝛽) ∑ 𝛽 𝑗 𝑦𝑡+𝑗+1
𝑗=0
𝑡
Thus, 𝜖𝑡 exactly reveals the gap between 𝑦𝑡 and 𝐸[𝑦𝑡 |𝑦+ ].
Next notice that representation (5.2), namely, 𝑦𝑡+1 − 𝑦𝑡 = −𝛽𝑎𝑡 + 𝑎𝑡+1 implies the linear difference equation
Solving this equation backward establishes that the one-step-prediction error 𝑎𝑡+1 is
∞
𝑎𝑡+1 = 𝑦𝑡+1 − (1 − 𝛽) ∑ 𝛽 𝑗 𝑦𝑡−𝑗 .
𝑗=0
Here the information set is 𝑦𝑡 = [𝑦𝑡 , 𝑦𝑡−1 , …] and a one step-ahead optimal prediction is
∞
𝐸[𝑦𝑡+1 |𝑦𝑡 ] = (1 − 𝛽) ∑ 𝛽 𝑗 𝑦𝑡−𝑗
𝑗=0
When we computed optimal consumption-saving policies for our two representations (5.1) and (5.2) by using formulas
obtained with the difference equation approach described in quantecon lecture, we obtained:
for a consumer having the information assumed in the news representation (5.1):
𝑐𝑡+1 − 𝑐𝑡 = 0
𝑏𝑡+1 − 𝑏𝑡 = −𝛽 −1 𝜖𝑡
for a consumer having the more limited information associated with the innovations representation (5.2):
𝑐𝑡+1 − 𝑐𝑡 = (1 − 𝛽 2 )𝑎𝑡+1
𝑏𝑡+1 − 𝑏𝑡 = −𝛽𝑎𝑡
These formulas agree with outcomes from Python programs below that deploy state-space representations and dynamic
programming.
Evidently, although they receive exactly the same histories of nonfinancial incomethe two consumers behave differently.
The better informed consumer who has the information sets associated with representation (5.1) responds to each shock
𝜖𝑡+1 by leaving his consumption unaltered and saving all of 𝜖𝑡+1 in anticipation of the permanently increased taxes that he
will bear in order to service the permanent interest payments on the risk-free bonds that the government has presumably
issued to pay for the one-time addition 𝜖𝑡+1 to his time 𝑡 + 1 nonfinancial income.
The less well informed consumer who has information sets associated with representation (5.2) responds to a shock 𝑎𝑡+1
by increasing his consumption by what he perceives to be the permanent part of the increase in consumption and by
increasing his saving by what he perceives to be the temporary part.
The behavior of the better informed consumer sharply illustrates the behavior predicted in a classic Ricardian equivalence
experiment.
We now cast our representations (5.1) and (5.2), respectively, in terms of the following two state space systems:
𝑦𝑡+1 1 −𝛽 −1 𝑦𝑡 𝜎
[ ]=[ ] [ ] + [ 𝜖 ] 𝑣𝑡+1
𝜖𝑡+1 0 0 𝜖𝑡 𝜎𝜖
(5.4)
𝑦
𝑦𝑡 = [1 0] [ 𝑡 ]
𝜖𝑡
and
𝑦𝑡+1 1 −𝛽 𝑦𝑡 𝜎
[ ]=[ ] [ ] + [ 𝑎 ] 𝑢𝑡+1
𝑎𝑡+1 0 0 𝑎𝑡 𝜎𝑎
(5.5)
𝑦
𝑦𝑡 = [1 0] [ 𝑡 ]
𝑎𝑡
where {𝑣𝑡 } and {𝑢𝑡 } are both i.i.d. sequences of univariate standardized normal random variables.
These two alternative income processes are ready to be used in the framework presented in the section “Comparison with
the Difference Equation Approach” in thid quantecon lecture.
All the code that we shall use below is presented in that lecture.
5.9 Computations
We shall use Python to form two state-space representations (5.4) and (5.5).
We set the following parameter values 𝜎𝜖 = 1, 𝜎𝑎 = 𝛽 −1 𝜎𝜖 = 𝛽 −1 where 𝛽 is the same value as the discount factor in
the household’s problem in the LQ savings problem in the lecture.
For these two representations, we use the code in this lecture to
• compute optimal decision rules for 𝑐𝑡 , 𝑏𝑡 for the two types of consumers associated with our two representations
of nonfinancial income
• use the value function objects 𝑃 , 𝑑 returned by the code to compute optimal values for the two representations
when evaluated at the initial condition
10
𝑥0 = [ ]
0
for each representation.
• create instances of the LinearStateSpace class for the two representations of the {𝑦𝑡 } process and use them to
obtain impulse response functions of 𝑐𝑡 and 𝑏𝑡 to the respective shocks 𝜖𝑡 and 𝑎𝑡 for the two representations.
• run simulations of {𝑦𝑡 , 𝑐𝑡 , 𝑏𝑡 } of length 𝑇 under both of the representations
We formulae the problem:
∞
2
min ∑ 𝛽 𝑡 (𝑐𝑡 − 𝛾)
𝑡=0
𝑦𝑡+1 1 −𝛽 −1 0 𝑦𝑡 0 𝜎𝜖
⎡ 𝜖 ⎤= ⎡ 0 0 0 ⎤ ⎡ 𝜖 ⎤ + ⎡ 0 ⎤ [ 𝑐 ] + ⎡ 𝜎 ⎤𝜈 ,
⎢ 𝑡+1 ⎥ ⎢ ⎥⎢ 𝑡 ⎥ ⎢ ⎥ 𝑡 ⎢ 𝜖 ⎥ 𝑡+1
⎣ − (1 + 𝑟)
⎣ 𝑏𝑡+1 ⎦ ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ 0 1 + 𝑟 ⎦ ⎣ 𝑏𝑡 ⎦ ⏟⎣⏟1⏟
+⏟𝑟⏟
⎦ ⎣
⏟ 0 ⎦
≡𝐴1 ≡𝐵1 ≡𝐶1
and
𝑦𝑡+1 1 −𝛽 0 𝑦𝑡 0 𝜎𝑎
⎡ 𝑎 ⎤ ⎡ 0 0 0 ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤
⎢ 𝑡+1 ⎥ = ⎢ ⎥ ⎢ 𝑎𝑡 ⎥ + ⎢ 0 ⎥ [ 𝑐𝑡 ] + ⎢ 𝜎𝑎 ⎥𝑢𝑡+1 .
⎣ 𝑏𝑡+1 ⎦ ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟
⎣ − (1 + 𝑟) 0 1 + 𝑟 ⎦ ⎣ 𝑏𝑡 ⎦ ⏟ ⎣⏟1⏟
+⏟𝑟⏟
⎦ ⎣ 0 ⎦
⏟
≡𝐴2 ≡𝐵2 ≡𝐶2
import numpy as np
import quantecon as qe
import matplotlib.pyplot as plt
5.9. Computations 83
Advanced Quantitative Economics with Python
# Set parameters
β, σϵ = 0.95, 1
σa = σϵ / β
R = 1 / β
Evidently, optimal consumption and debt decision rules for the consumer having news representation (5.1) are
𝑐𝑡∗ = 𝑦𝑡 − 𝜖𝑡 − (1 − 𝛽) 𝑏𝑡 ,
∗
𝑏𝑡+1 = 𝛽 −1 𝑐𝑡∗ + 𝛽 −1 𝑏𝑡 − 𝛽 −1 𝑦𝑡
= 𝛽 −1 𝑦𝑡 − 𝛽 −1 𝜖𝑡 − (𝛽 −1 − 1) 𝑏𝑡 + 𝛽 −1 𝑏𝑡 − 𝛽 −1 𝑦𝑡
= 𝑏𝑡 − 𝛽 −1 𝜖𝑡 .
# Innovations representation
ALQ2 = np.array([[1, -β, 0],
[0, 0, 0],
[-R, 0, R]])
BLQ2 = np.array([[0, 0, R]]).T
CLQ2 = np.array([[σa, σa, 0]]).T
-F2
For a consumer having access only to the information associated with the innovations representation (5.2), the optimal
Now we construct two Linear State Space models that emerge from using optimal policies of the form 𝑢𝑡 = −𝐹 𝑥𝑡 .
Take the more informative original representation (5.1) as an example:
𝑦𝑡+1 𝑦𝑡
⎡ 𝜖 ⎤ = (𝐴 − 𝐵 𝐹 ) ⎡ 𝜖 ⎤ + 𝐶 𝜈
⎢ 𝑡+1 ⎥ 1 1 1 ⎢ 𝑡 ⎥ 1 𝑡+1
⎣ 𝑏𝑡+1 ⎦ ⎣ 𝑏𝑡 ⎦
𝑦
𝑐𝑡 −𝐹1 ⎡ 𝑡 ⎤
[ ]=[ ] ⎢ 𝜖𝑡 ⎥
𝑏𝑡 𝑆𝑏
⎣ 𝑏𝑡 ⎦
To have the Linear State Space model be of an innovations representation form (5.2), we can simply replace the corre-
sponding matrices.
5.9. Computations 85
Advanced Quantitative Economics with Python
<matplotlib.legend.Legend at 0x7fc14e8fe450>
The above two impulse response functions show that when the consumer has the information assumed in the more infor-
mative representation (5.1), his response to receiving a positive shock of 𝜖𝑡 is to leave his consumption unchanged and to
save the entire amount of his extra income and then forever roll over the extra bonds that he holds.
To see this notice, that starting from next period on, his debt permanently decreases by 𝛽 −1
plt.title("innovations representation")
plt.plot(range(J), c_res2 / σa, label="c impulse response function")
plt.plot(range(J), b_res2 / σa, label="b impulse response function")
plt.plot([0, J-1], [0, 0], '--', color='k')
plt.legend()
<matplotlib.legend.Legend at 0x7fc14e6a64e0>
The above impulse responses show that when the consumer has only the information that is assumed to be available
under the innovations representation (5.2) for {𝑦𝑡 − 𝑦𝑡−1 }, he responds to a positive 𝑎𝑡 by permanently increasing his
consumption.
He accomplishes this by consuming a fraction (1 − 𝛽 2 ) of the increment 𝑎𝑡 to his nonfinancial income and saving the
rest, thereby lowering 𝑏𝑡+1 in order to finance the permanent increment in his consumption.
The preceding computations confirm what we had derived earlier using paper and pencil.
Now let’s simulate some paths of consumption and debt for our two types of consumers while always presenting both
types with the same {𝑦𝑡 } path.
x1, y1 = LSS1.simulate(ts_length=T)
plt.plot(range(T), y1[0, :], label="c")
plt.plot(range(T), x1[2, :], label="b")
plt.plot(range(T), x1[0, :], label="y")
plt.title("more informative representation")
plt.legend()
<matplotlib.legend.Legend at 0x7fc14e6a6d20>
5.9. Computations 87
Advanced Quantitative Economics with Python
x2, y2 = LSS2.simulate(ts_length=T)
plt.plot(range(T), y2[0, :], label="c")
plt.plot(range(T), x2[2, :], label="b")
plt.plot(range(T), x2[0, :], label="y")
plt.title("innovations representation")
plt.legend()
<matplotlib.legend.Legend at 0x7fc14e1b6c00>
We now form a single {𝑦𝑡 }𝑇𝑡=0 realization that we will use to simulate decisions associated with our two types of consumer.
We accomplish this in the following steps.
1. We form a {𝑦𝑡 , 𝜖𝑡 } realization by drawing a long simulation of {𝜖𝑡 }𝑇𝑡=0 , where 𝑇 is a big integer, 𝜖𝑡 = 𝜎𝜖 𝑣𝑡 , 𝑣𝑡 is
a standard normal scalar, 𝑦0 = 100, and
𝑦𝑡+1 − 𝑦𝑡 = −𝛽 −1 𝜖𝑡 + 𝜖𝑡+1 .
2. We take the {𝑦𝑡 } realization generated in step 1 and form an innovation process {𝑎𝑡 } from the formulas
𝑎0 = 0
𝑡−1
𝑎𝑡 = ∑ 𝛽 𝑗 (𝑦𝑡−𝑗 − 𝑦𝑡−𝑗−1 ) + 𝛽 𝑡 𝑎0 , 𝑡≥1
𝑗=0
3. We throw away the first 𝑆 observations and form a sample {𝑦𝑡 , 𝜖𝑡 , 𝑎𝑡 }𝑇𝑆+1 as the realization that we’ll use in the
following steps.
4. We use the step 3 realization to evaluate and simulate the decision rules for 𝑐𝑡 , 𝑏𝑡 that Python has computed for
us above.
The above steps implement the experiment of comparing decisions made by two consumers having identical incomes at
each date but at each date having different information about their future incomes.
Here we use formula (5.3) above to compute 𝑎𝑡+1 as a function of the history 𝜖𝑡+1 , 𝜖𝑡 , 𝜖𝑡−1 , …
Thus, we compute
We can verify that we recover the same {𝑎𝑡 } sequence computed earlier.
This quantecon lecture contains another example of a shock-invertibility issue that is endemic to the LQ permanent income
or consumption smoothing model.
The technical issue discussed there is ultimately the source of the shock-invertibility issues discussed by Eric Leeper,
Todd Walker, and Susan Yang [Leeper et al., 2013] in their analysis of fiscal foresight.
SIX
6.1 Overview
91
Advanced Quantitative Economics with Python
import numpy as np
import quantecon as qe
import matplotlib.pyplot as plt
import scipy.linalg as la
6.2 Background
We’ll study a complete markets model adapted to a setting with a continuous Markov state like that in the first lecture on
the permanent income model.
In that model
• a consumer can trade only a single risk-free one-period bond bearing gross one-period risk-free interest rate equal
to 𝛽 −1 .
• a consumer’s exogenous nonfinancial income is governed by a linear state space model driven by Gaussian shocks,
the kind of model studied in an earlier lecture about linear state space models.
Let’s write down a complete markets counterpart of that model.
Suppose that nonfinancial income is governed by the state space system
where 𝜙(⋅ | 𝜇, Σ) is a multivariate Gaussian distribution with mean vector 𝜇 and covariance matrix Σ.
With the pricing kernel 𝑞𝑡+1 (𝑥𝑡+1 | 𝑥𝑡 ) in hand, we can price claims to consumption at time 𝑡 + 1 consumption that pay
off when 𝑥𝑡+1 ∈ 𝑆 at time 𝑡 + 1:
𝑛
where 𝑆 is a subset of ℝ .
The price ∫𝑆 𝑞𝑡+1 (𝑥𝑡+1 | 𝑥𝑡 )𝑑𝑥𝑡+1 of such a claim depends on state 𝑥𝑡 because the prices of the 𝑥𝑡+1 -contingent securities
depend on 𝑥𝑡 through the pricing kernel 𝑞(𝑥𝑡+1 | 𝑥𝑡 ).
Let 𝑏(𝑥𝑡+1 ) be a vector of state-contingent debt due at 𝑡 + 1 as a function of the 𝑡 + 1 state 𝑥𝑡+1 .
Using the pricing kernel assumed in (6.1), the value at 𝑡 of 𝑏(𝑥𝑡+1 ) is evidently
In our complete markets setting, the consumer faces a sequence of budget constraints
or
which verifies that 𝛽𝐸𝑡 𝑏𝑡+1 is the value of time 𝑡 + 1 state-contingent claims on time 𝑡 + 1 consumption issued by the
consumer at time 𝑡
We can solve the time 𝑡 budget constraint forward to obtain
∞
𝑏𝑡 = 𝔼𝑡 ∑ 𝛽 𝑗 (𝑦𝑡+𝑗 − 𝑐𝑡+𝑗 )
𝑗=0
In the incomplete markets version of the model, we assumed that 𝑢(𝑐𝑡 ) = −(𝑐𝑡 − 𝛾)2 , so that the above utility functional
became
∞
− ∑ 𝛽 𝑡 (𝑐𝑡 − 𝛾)2 , 0<𝛽<1
𝑡=0
But in the complete markets version, it is tractable to assume a more general utility function that satisfies 𝑢′ > 0 and
𝑢″ < 0.
First-order conditions for the consumer’s problem with complete markets and our assumption about Arrow securities
prices are
or
1
𝑏𝑡 = 𝑆𝑦 (𝐼 − 𝛽𝐴)−1 𝑥𝑡 − 𝑐̄ (6.2)
1−𝛽
where 𝑐 ̄ satisfies
1
𝑏̄0 = 𝑆𝑦 (𝐼 − 𝛽𝐴)−1 𝑥0 − 𝑐̄ (6.3)
1−𝛽
where 𝑏̄0 is an initial level of the consumer’s debt due at time 𝑡 = 0, specified as a parameter of the problem.
Thus, in the complete markets version of the consumption-smoothing model, 𝑐𝑡 = 𝑐,̄ ∀𝑡 ≥ 0 is determined by (6.3) and
the consumer’s debt is the fixed function of the state 𝑥𝑡 described by (6.2).
Please recall that in the LQ permanent income model studied in permanent income model, the state is 𝑥𝑡 , 𝑏𝑡 , where 𝑏𝑡 is
a complicated function of past state vectors 𝑥𝑡−𝑗 .
Notice that in contrast to that incomplete markets model, at time 𝑡 the state vector is 𝑥𝑡 alone in our complete markets
model.
Here’s an example that shows how in this setting the availability of insurance against fluctuating nonfinancial income
allows the consumer completely to smooth consumption across time and across states of the world
# Debt
x_hist, y_hist = lss.simulate(T)
b_hist = np.squeeze(S_y @ rm @ x_hist - cbar / (1 - β))
# Define parameters
N_simul = 80
α, ρ1, ρ2 = 10.0, 0.9, 0.0
σ = 1.0
# Consumption plots
ax[0].set_title('Consumption and income')
ax[0].plot(np.arange(N_simul), c_hist_com, label='consumption')
ax[0].plot(np.arange(N_simul), y_hist_com, label='income', alpha=.6, linestyle='--')
ax[0].legend()
ax[0].set_xlabel('Periods')
ax[0].set_ylim([80, 120])
# Debt plots
ax[1].set_title('Debt and income')
ax[1].plot(np.arange(N_simul), b_hist_com, label='debt')
ax[1].plot(np.arange(N_simul), y_hist_com, label='Income', alpha=.6, linestyle='--')
ax[1].legend()
ax[1].axhline(0, color='k')
ax[1].set_xlabel('Periods')
plt.show()
The incomplete markets version of the model with nonfinancial income being governed by a linear state space system is
described in permanent income model.
In that incomplete markerts setting, consumption follows a random walk and the consumer’s debt follows a process with
a unit root.
We now turn to a finite-state Markov version of the model in which the consumer’s nonfinancial income is an exact
function of a Markov state that takes one of 𝑁 values.
We’ll start with a setting in which in each version of our consumption-smoothing model, nonfinancial income is governed
by a two-state Markov chain (it’s easy to generalize this to an 𝑁 state Markov chain).
In particular, the state 𝑠𝑡 ∈ {1, 2} follows a Markov chain with transition probability matrix
𝑃𝑖𝑗 = ℙ{𝑠𝑡+1 = 𝑗 | 𝑠𝑡 = 𝑖}
𝑦1̄ if 𝑠𝑡 = 1
𝑦𝑡 = {
𝑦2̄ if 𝑠𝑡 = 2
Our complete and incomplete markets models differ in how thoroughly the market structure allows a consumer to transfer
resources across time and Markov states, there being more transfer opportunities in the complete markets setting than in
the incomplete markets setting.
Watch how these differences in opportunities affect
• how smooth consumption is across time and Markov states
• how the consumer chooses to make his levels of indebtedness behave over time and across Markov states
At each date 𝑡 ≥ 0, the consumer trades a full array of one-period ahead Arrow securities.
We assume that prices of these securities are exogenous to the consumer.
Exogenous means that they are unaffected by the consumer’s decisions.
In Markov state 𝑠𝑡 at time 𝑡, one unit of consumption in state 𝑠𝑡+1 at time 𝑡 + 1 costs 𝑞(𝑠𝑡+1 | 𝑠𝑡 ) units of the time 𝑡
consumption good.
The prices 𝑞(𝑠𝑡+1 | 𝑠𝑡 ) are given and can be organized into a matrix 𝑄 with 𝑄𝑖𝑗 = 𝑞(𝑗|𝑖)
At time 𝑡 = 0, the consumer starts with an inherited level of debt due at time 0 of 𝑏0 units of time 0 consumption goods.
The consumer’s budget constraint at 𝑡 ≥ 0 in Markov state 𝑠𝑡 is
where 𝑏𝑡 is the consumer’s one-period debt that falls due at time 𝑡 and 𝑏𝑡+1 (𝑗 | 𝑠𝑡 ) are the consumer’s time 𝑡 sales of the
time 𝑡 + 1 consumption good in Markov state 𝑗.
Thus
• 𝑞(𝑗 | 𝑠𝑡 )𝑏𝑡+1 (𝑗 | 𝑠𝑡 ) is a source of time 𝑡 financial income for the consumer in Markov state 𝑠𝑡
• 𝑏𝑡 ≡ 𝑏𝑡 (𝑗 | 𝑠𝑡−1 ) is a source of time 𝑡 expenditures for the consumer when 𝑠𝑡 = 𝑗
Remark: We are ignoring an important technicality here, namely, that the consumer’s choice of 𝑏𝑡+1 (𝑗| 𝑠𝑡 ) must respect
so-called natural debt limits that assure that it is feasible for the consumer to repay debts due even if he consumers zero
forevermore. We shall discuss such debt limits in another lecture.
A natural analog of Hall’s assumption that the one-period risk-free gross interest rate is 𝛽 −1 is
To understand how this is a natural analogue, observe that in state 𝑖 it costs ∑𝑗 𝑞(𝑗 | 𝑖) to purchase one unit of consumption
next period for sure, i.e., meaning no matter what Markov state 𝑗 occurs at 𝑡 + 1.
Hence the implied price of a risk-free claim on one unit of consumption next period is
∑ 𝑞(𝑗 | 𝑖) = ∑ 𝛽𝑃𝑖𝑗 = 𝛽
𝑗 𝑗
This confirms the sense in which (6.6) is a natural counterpart to Hall’s assumption that the risk-free one-period gross
interest rate is 𝑅 = 𝛽 −1 .
It is timely please to recall that the gross one-period risk-free interest rate is the reciprocal of the price at time 𝑡 of a
risk-free claim on one unit of consumption tomorrow.
First-order necessary conditions for maximizing the consumer’s expected utility subject to the sequence of budget con-
straints (6.5) are
𝑢′ (𝑐𝑡+1 )
𝛽 ℙ{𝑠𝑡+1 | 𝑠𝑡 } = 𝑞(𝑠𝑡+1 | 𝑠𝑡 )
𝑢′ (𝑐𝑡 )
for all 𝑠𝑡 , 𝑠𝑡+1 or, under our assumption (6.6) about Arrow security prices,
𝑐𝑡+1 = 𝑐𝑡 (6.7)
Thus, our consumer sets 𝑐𝑡 = 𝑐 ̄ for all 𝑡 ≥ 0 for some value 𝑐 ̄ that it is our job now to determine along with values for
𝑏𝑡+1 (𝑗|𝑠𝑡 = 𝑖) for 𝑖 = 1, 2 and 𝑗 = 1, 2.
We’ll use a guess and verify method to determine these objects
Guess: We’ll make the plausible guess that
so that the amount borrowed today depends only on tomorrow’s Markov state. (Why is this is a plausible guess?)
To determine 𝑐,̄ we shall deduce implications of the consumer’s budget constraints in each Markov state today and our
guess (6.8) about the consumer’s debt level choices.
For 𝑡 ≥ 1, these imply
where 𝑏0 is the (exogenous) debt the consumer is assumed to bring into period 0
If we substitute (6.10) into the first equation of (6.9) and rearrange, we discover that
𝑏(1) = 𝑏0 (6.11)
We can then use the second equation of (6.9) to deduce the restriction
The preceding calculations indicate that in the complete markets version of our model, we obtain the following striking
results:
• The consumer chooses to make consumption perfectly constant across time and across Markov states.
• State-contingent debt purchases 𝑏𝑡+1 (𝑠𝑡+1 = 𝑗|𝑠𝑡 = 𝑖) depend only on 𝑗
• If the initial Markov state is 𝑠0 = 𝑗 and initial consumer debt is 𝑏0 , then debt in Markov state 𝑗 satisfies 𝑏(𝑗) = 𝑏0
To summarize what we have achieved up to now, we have computed the constant level of consumption 𝑐 ̄ and indicated
how that level depends on the underlying specifications of preferences, Arrow securities prices, the stochastic process of
exogenous nonfinancial income, and the initial debt level 𝑏0
• The consumer’s debt neither accumulates, nor decumulates, nor drifts – instead, the debt level each period is an
exact function of the Markov state, so in the two-state Markov case, it switches between two values.
• We have verified guess (6.8).
• When the state 𝑠𝑡 returns to the initial state 𝑠0 , debt returns to the initial debt level.
• Debt levels in all other states depend on virtually all remaining parameters of the model.
6.4.2 Code
Here’s some code that, among other things, contains a function called consumption_complete().
This function computes {𝑏(𝑖)}𝑁
𝑖=1 , 𝑐 ̄ as outcomes given a set of parameters for the general case with 𝑁 Markov states
under the assumption of complete markets
class ConsumptionProblem:
"""
The data for a consumption problem, including some default values.
"""
def __init__(self,
β=.96,
y=[2, 1.5],
b0=3,
P=[[.8, .2],
[.4, .6]],
init=0):
"""
Parameters
----------
β : discount factor
y : list containing the two income levels
b0 : debt in period 0 (= initial state debt level)
(continues on next page)
return s_path
def consumption_complete(cp):
"""
Computes endogenous values for the complete market case.
Parameters
----------
cp : instance of ConsumptionProblem
Returns
-------
Q = β * P
"""
β, P, y, b0, init = cp.β, cp.P, cp.y, cp.b0, cp.init # Unpack
A = np.zeros((n, n))
A[:, 0] = 1
A[1:, 1:] = np.eye(n-1)
c_bar = x[0, 0]
b = x[1:, 0]
return c_bar, b
Parameters
----------
cp : instance of ConsumptionProblem
s_path : the path of states
"""
β, P, y, b0 = cp.β, cp.P, cp.y, cp.b0 # Unpack
N_simul = len(s_path)
# Useful variables
n = len(y)
y.shape = (n, 1)
v = np.linalg.inv(np.eye(n) - β * P) @ y
for i, s in enumerate(s_path):
c_path[i] = (1 - β) * (v - np.full((n, 1), b_path[i]))[s, 0]
b_path[i + 1] = b_path[i] + db[s, 0]
cp = ConsumptionProblem()
c_bar, b = consumption_complete(cp)
np.isclose(c_bar + b[1] - cp.y[1] - (cp.β * cp.P)[1, :] @ b, 0)
True
Below, we’ll take the outcomes produced by this code – in particular the implied consumption and debt paths – and
compare them with outcomes from an incomplete markets model in the spirit of Hall [Hall, 1978]
This is a version of the original model of Hall (1978) in which the consumer’s ability to substitute intertemporally is
constrained by his ability to buy or sell only one security, a risk-free one-period bond bearing a constant gross interest
rate that equals 𝛽 −1 .
Given an initial debt 𝑏0 at time 0, the consumer faces a sequence of budget constraints
𝑐𝑡 + 𝑏𝑡 = 𝑦𝑡 + 𝛽𝑏𝑡+1 , 𝑡≥0
where 𝛽 is the price at time 𝑡 of a risk-free claim on one unit of time consumption at time 𝑡 + 1.
First-order conditions for the consumer’s problem are
which for our finite-state Markov setting is Hall’s (1978) conclusion that consumption follows a random walk.
As we saw in our first lecture on the permanent income model, this leads to
∞
𝑏𝑡 = 𝔼𝑡 ∑ 𝛽 𝑗 𝑦𝑡+𝑗 − (1 − 𝛽)−1 𝑐𝑡 (6.14)
𝑗=0
and
∞
𝑐𝑡 = (1 − 𝛽) [𝔼𝑡 ∑ 𝛽 𝑗 𝑦𝑡+𝑗 − 𝑏𝑡 ] (6.15)
𝑗=0
Equation (6.15) expresses 𝑐𝑡 as a net interest rate factor 1 − 𝛽 times the sum of the expected present value of nonfinancial
∞
income 𝔼𝑡 ∑𝑗=0 𝛽 𝑗 𝑦𝑡+𝑗 and financial wealth −𝑏𝑡 .
Substituting (6.15) into the one-period budget constraint and rearranging leads to
∞
𝑏𝑡+1 − 𝑏𝑡 = 𝛽 −1 [(1 − 𝛽)𝔼𝑡 ∑ 𝛽 𝑗 𝑦𝑡+𝑗 − 𝑦𝑡 ] (6.16)
𝑗=0
∞
Now let’s calculate the key term 𝔼𝑡 ∑𝑗=0 𝛽 𝑗 𝑦𝑡+𝑗 in our finite Markov chain setting.
Define the expected discounted present value of non-financial income
∞
𝑣𝑡 ∶= 𝔼𝑡 ∑ 𝛽 𝑗 𝑦𝑡+𝑗
𝑗=0
𝑣𝑡 ∶= 𝑦𝑡 + 𝛽𝔼𝑡 𝑣𝑡+1
In our two-state Markov chain setting, 𝑣𝑡 = 𝑣(1) when 𝑠𝑡 = 1 and 𝑣𝑡 = 𝑣(2) when 𝑠𝑡 = 2.
Therefore, we can write our Bellman equation as
or
𝑣 ⃗ = 𝑦 ⃗ + 𝛽𝑃 𝑣 ⃗
𝑣(1) 𝑦(1)
where 𝑣 ⃗ = [ ] and 𝑦 ⃗ = [ ].
𝑣(2) 𝑦(2)
We can also write the last expression as
𝑣 ⃗ = (𝐼 − 𝛽𝑃 )−1 𝑦 ⃗
In our finite Markov chain setting, from expression (6.15), consumption at date 𝑡 when debt is 𝑏𝑡 and the Markov state
today is 𝑠𝑡 = 𝑖 is evidently
In contrast to outcomes in the complete markets model, in the incomplete markets model
• consumption drifts over time as a random walk; the level of consumption at time 𝑡 depends on the level of debt that
the consumer brings into the period as well as the expected discounted present value of nonfinancial income at 𝑡.
• the consumer’s debt drifts upward over time in response to low realizations of nonfinancial income and drifts
downward over time in response to high realizations of nonfinancial income.
• the drift over time in the consumer’s debt and the dependence of current consumption on today’s debt level account
for the drift over time in consumption.
The code above also contains a function called consumption_incomplete() that uses (6.17) and (6.18) to
• simulate paths of 𝑦𝑡 , 𝑐𝑡 , 𝑏𝑡+1
• plot these against values of 𝑐,̄ 𝑏(𝑠1 ), 𝑏(𝑠2 ) found in a corresponding complete markets economy
Let’s try this, using the same parameters in both complete and incomplete markets economies
cp = ConsumptionProblem()
s_path = cp.simulate()
N_simul = len(s_path)
ax[0].set_title('Consumption paths')
ax[0].plot(np.arange(N_simul), c_path, label='incomplete market')
ax[0].plot(np.arange(N_simul), np.full(N_simul, c_bar),
(continues on next page)
ax[1].set_title('Debt paths')
ax[1].plot(np.arange(N_simul), debt_path, label='incomplete market')
ax[1].plot(np.arange(N_simul), debt_complete[s_path],
label='complete market')
ax[1].plot(np.arange(N_simul), y_path, label='income', alpha=.6, ls='--')
ax[1].legend()
ax[1].axhline(0, color='k', ls='--')
ax[1].set_xlabel('Periods')
plt.show()
In the graph on the left, for the same sample path of nonfinancial income 𝑦𝑡 , notice that
• consumption is constant when there are complete markets, but takes a random walk in the incomplete markets
version of the model.
• the consumer’s debt oscillates between two values that are functions of the Markov state in the complete markets
model, while the consumer’s debt drifts in a “unit root” fashion in the incomplete markets economy.
6.5.3 A sequel
In tax smoothing with complete and incomplete markets, we reinterpret the mathematics and Python code presented in this
lecture in order to construct tax-smoothing models in the incomplete markets tradition of Barro [Barro, 1979] as well as
in the complete markets tradition of Lucas and Stokey [Lucas and Stokey, 1983].
SEVEN
7.1 Overview
This lecture describes tax-smoothing models that are counterparts to consumption-smoothing models in Consumption
Smoothing with Complete and Incomplete Markets.
• one is in the complete markets tradition of Lucas and Stokey [Lucas and Stokey, 1983].
• the other is in the incomplete markets tradition of Barro [Barro, 1979].
Complete markets allow a government to buy or sell claims contingent on all possible Markov states.
Incomplete markets allow a government to buy or sell only a limited set of securities, often only a single risk-free security.
Barro [Barro, 1979] worked in an incomplete markets tradition by assuming that the only asset that can be traded is a
risk-free one period bond.
In his consumption-smoothing model, Hall [Hall, 1978] had assumed an exogenous stochastic process of nonfinancial
income and an exogenous gross interest rate on one period risk-free debt that equals 𝛽 −1 , where 𝛽 ∈ (0, 1) is also a
consumer’s intertemporal discount factor.
Barro [Barro, 1979] made an analogous assumption about the risk-free interest rate in a tax-smoothing model that turns
out to have the same mathematical structure as Hall’s consumption-smoothing model.
To get Barro’s model from Hall’s, all we have to do is to rename variables.
We maintain Hall’s and Barro’s assumption about the interest rate when we describe an incomplete markets version of
our model.
In addition, we extend their assumption about the interest rate to an appropriate counterpart to create a “complete markets”
model in the style of Lucas and Stokey [Lucas and Stokey, 1983].
105
Advanced Quantitative Economics with Python
For each version of a consumption-smoothing model, a tax-smoothing counterpart can be obtained simply by relabeling
• consumption as tax collections
• a consumer’s one-period utility function as a government’s one-period loss function from collecting taxes that im-
pose deadweight welfare losses
• a consumer’s nonfinancial income as a government’s purchases
• a consumer’s debt as a government’s assets
Thus, we can convert the consumption-smoothing models in lecture Consumption Smoothing with Complete and Incomplete
Markets into tax-smoothing models by setting 𝑐𝑡 = 𝑇𝑡 , 𝑦𝑡 = 𝐺𝑡 , and −𝑏𝑡 = 𝑎𝑡 , where 𝑇𝑡 is total tax collections, {𝐺𝑡 }
is an exogenous government expenditures process, and 𝑎𝑡 is the government’s holdings of one-period risk-free bonds
coming maturing at the due at the beginning of time 𝑡.
For elaborations on this theme, please see Optimal Savings II: LQ Techniques and later parts of this lecture.
We’ll spend most of this lecture studying acquire finite-state Markov specification, but will also treat the linear state space
specification.
Link to History
For those who love history, President Thomas Jefferson’s Secretary of Treasury Albert Gallatin (1807) [Gallatin, 1837]
seems to have prescribed policies that come from Barro’s model [Barro, 1979]
Let’s start with some standard imports:
import numpy as np
import quantecon as qe
import matplotlib.pyplot as plt
To exploit the isomorphism between consumption-smoothing and tax-smoothing models, we simply use code from Con-
sumption Smoothing with Complete and Incomplete Markets
7.1.2 Code
class ConsumptionProblem:
"""
The data for a consumption problem, including some default values.
"""
def __init__(self,
β=.96,
y=[2, 1.5],
b0=3,
P=[[.8, .2],
[.4, .6]],
init=0):
"""
(continues on next page)
β : discount factor
y : list containing the two income levels
b0 : debt in period 0 (= initial state debt level)
P : 2x2 transition matrix
init : index of initial state s0
"""
self.β = β
self.y = np.asarray(y)
self.b0 = b0
self.P = np.asarray(P)
self.init = init
return s_path
def consumption_complete(cp):
"""
Computes endogenous values for the complete market case.
Parameters
----------
cp : instance of ConsumptionProblem
Returns
-------
Q = β * P
"""
β, P, y, b0, init = cp.β, cp.P, cp.y, cp.b0, cp.init # Unpack
A = np.zeros((n, n))
A[:, 0] = 1
A[1:, 1:] = np.eye(n-1)
c_bar = x[0, 0]
b = x[1:, 0]
return c_bar, b
Parameters
----------
cp : instance of ConsumptionProblem
s_path : the path of states
"""
β, P, y, b0 = cp.β, cp.P, cp.y, cp.b0 # Unpack
N_simul = len(s_path)
# Useful variables
n = len(y)
y.shape = (n, 1)
v = np.linalg.inv(np.eye(n) - β * P) @ y
for i, s in enumerate(s_path):
c_path[i] = (1 - β) * (v - np.full((n, 1), b_path[i]))[s, 0]
b_path[i + 1] = b_path[i] + db[s, 0]
The code above also contains a function called consumption_incomplete() that uses (6.17) and (6.18) to
• simulate paths of 𝑦𝑡 , 𝑐𝑡 , 𝑏𝑡+1
• plot these against values of 𝑐,̄ 𝑏(𝑠1 ), 𝑏(𝑠2 ) found in a corresponding complete markets economy
Let’s try this, using the same parameters in both complete and incomplete markets economies
cp = ConsumptionProblem()
s_path = cp.simulate()
N_simul = len(s_path)
ax[0].set_title('Consumption paths')
ax[0].plot(np.arange(N_simul), c_path, label='incomplete market')
ax[0].plot(np.arange(N_simul), np.full(N_simul, c_bar), label='complete market')
ax[0].plot(np.arange(N_simul), y_path, label='income', alpha=.6, ls='--')
ax[0].legend()
ax[0].set_xlabel('Periods')
ax[1].set_title('Debt paths')
ax[1].plot(np.arange(N_simul), debt_path, label='incomplete market')
ax[1].plot(np.arange(N_simul), debt_complete[s_path], label='complete market')
ax[1].plot(np.arange(N_simul), y_path, label='income', alpha=.6, ls='--')
ax[1].legend()
ax[1].axhline(0, color='k', ls='--')
ax[1].set_xlabel('Periods')
plt.show()
In the graph on the left, for the same sample path of nonfinancial income 𝑦𝑡 , notice that
• consumption is constant when there are complete markets.
• consumption takes a random walk in the incomplete markets version of the model.
• the consumer’s debt oscillates between two values that are functions of the Markov state in the complete markets
model.
• the consumer’s debt drifts because it contains a unit root in the incomplete markets economy.
As indicated above, we relabel variables to acquire tax-smoothing interpretations of the complete markets and incomplete
markets consumption-smoothing models.
plt.show()
𝑇𝑖 + 𝑏𝑖 = 𝐺𝑖 + ∑ 𝑄𝑖𝑗 𝑏𝑗
𝑗
where
𝑄𝑖𝑗 = 𝛽𝑃𝑖𝑗
is the price today of one unit of goods in Markov state 𝑗 tomorrow when the Markov state is 𝑖 today.
𝑏𝑖 is the government’s level of assets when it arrives in Markov state 𝑖.
That is, 𝑏𝑖 equals one-period state-contingent claims owed to the government that fall due at time 𝑡 when the Markov state
is 𝑖.
Thus, if 𝑏𝑖 < 0, it means the government is owed 𝑏𝑖 or owes −𝑏𝑖 when the economy arrives in Markov state 𝑖 at time 𝑡.
In our examples below, this happens when in a previous war-time period the government has sold an Arrow securities
paying off −𝑏𝑖 in peacetime Markov state 𝑖
It can be enlightening to express the government’s budget constraint in Markov state 𝑖 as
𝑇𝑖 = 𝐺𝑖 + (∑ 𝑄𝑖𝑗 𝑏𝑗 − 𝑏𝑖 )
𝑗
in which the term (∑𝑗 𝑄𝑖𝑗 𝑏𝑗 − 𝑏𝑖 ) equals the net amount that the government spends to purchase one-period Arrow
securities that will pay off next period in Markov states 𝑗 = 1, … , 𝑁 after it has received payments 𝑏𝑖 this period.
The cumulative return earned from putting 1 unit of time 𝑡 goods into the government portfolio of state-contingent
securities at time 𝑡 and then rolling over the proceeds into the government portfolio each period thereafter is
Here is some code that computes one-period and cumulative returns on the government portfolio in the finite-state Markov
version of our complete markets model.
Convention: In this code, when 𝑃𝑖𝑗 = 0, we arbitrarily set 𝑅(𝑗|𝑖) to be 0.
values = Q @ b
n = len(b)
R = np.zeros((n, n))
for i in range(n):
ind = cp.P[i, :] != 0
R[i, ind] = b[ind] / values[i]
return R
RT_path = np.empty(T)
RT_path[0] = 1
RT_path[1:] = np.cumprod([R[s_path[t], s_path[t+1]] for t in range(T-1)])
return RT_path
# Parameters
β = .96
cp = ConsumptionProblem(β, g, b0, P)
(continues on next page)
print(f"P \n {P}")
print(f"Q \n {Q}")
print(f"Govt expenditures in peace and war = {g}")
print(f"Constant tax collections = {T_bar}")
print(f"Govt debts in two states = {-b}")
msg = """
Now let's check the government's budget constraint in peace and war.
Our assumptions imply that the government always purchases 0 units of the
Arrow peace security.
"""
print(msg)
AS1 = Q[0, :] @ b
# spending on Arrow security
# since the spending on Arrow peace security is not 0 anymore after we change b0 to 1
print(f"Spending on Arrow security in peace = {AS1}")
AS2 = Q[1, :] @ b
print(f"Spending on Arrow security in war = {AS2}")
print("")
# tax collections minus debt levels
print("Government tax collections minus debt levels in peace and war")
TB1 = T_bar + b[0]
print(f"T+b in peace = {TB1}")
TB2 = T_bar + b[1]
print(f"T+b in war = {TB2}")
print("")
print("Total government spending in peace and war")
G1 = g[0] + AS1
G2 = g[1] + AS2
print(f"Peace = {G1}")
print(f"War = {G2}")
print("")
print("Let's see ex-post and ex-ante returns on Arrow securities")
Π = np.reciprocal(Q)
exret = Π
print(f"Ex-post returns to purchase of Arrow securities = \n {exret}")
exant = Π * P
print(f"Ex-ante returns to purchase of Arrow securities \n {exant}")
print("")
print("The Ex-post one-period gross return on the portfolio of government assets")
print(R)
print(RT_path[-1])
P
[[0.8 0.2]
[0.4 0.6]]
Q
[[0.768 0.192]
[0.384 0.576]]
Govt expenditures in peace and war = [1, 2]
Constant tax collections = 1.2716883116883118
Govt debts in two states = [-1. -2.62337662]
Now let's check the government's budget constraint in peace and war.
Our assumptions imply that the government always purchases 0 units of the
Arrow peace security.
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
2.0860704239993675
7.3.2 Explanation
In this example, the government always purchase 1 units of the Arrow security that pays off in peace time (Markov state
1).
And it purchases a higher amount of the security that pays off in war time (Markov state 2).
Thus, this is an example in which
• during peacetime, the government purchases insurance against the possibility that war breaks out next period
• during wartime, the government purchases insurance against the possibility that war continues another period
• so long as peace continues, the ex post return on insurance against war is low
• when war breaks out or continues, the ex post return on insurance against war is high
• given the history of states that we assumed, the value of one unit of the portfolio of government assets eventually
doubles in the end because of high returns during wartime.
We recommend plugging the quantities computed above into the government budget constraints in the two Markov states
and staring.
Exercise 7.3.1
Try changing the Markov transition matrix so that
1 0
𝑃 =[ ]
.2 .8
Also, start the system in Markov state 2 (war) with initial government assets −10, so that the government starts the war
in debt and 𝑏2 = −10.
To interpret some episodes in the fiscal history of the United States, we find it interesting to study a few more examples.
We compute examples in an 𝑁 state Markov setting under both complete and incomplete markets.
These examples differ in how Markov states are jumping between peace and war.
To wrap procedures for solving models, relabeling graphs so that we record government debt rather than government
assets, and displaying results, we construct a Python class.
class TaxSmoothingExample:
"""
construct a tax-smoothing example, by relabeling consumption problem class.
"""
def __init__(self, g, P, b0, states, β=.96,
init=0, s_path=None, N_simul=80, random_state=1):
def display(self):
# plot graphs
N = len(self.T_path)
plt.figure()
plt.title('Tax collection paths')
plt.plot(np.arange(N), self.T_path, label='incomplete market')
plt.plot(np.arange(N), np.full(N, self.T_bar), label='complete market')
plt.plot(np.arange(N), self.g_path, label='govt expenditures', alpha=.6, ls='-
↪-')
plt.legend()
plt.xlabel('Periods')
plt.show()
plt.legend()
plt.axhline(0, color='k', ls='--')
plt.xlabel('Periods')
plt.show()
fig, ax = plt.subplots()
ax.set_title('Cumulative return path (complete markets)')
line1 = ax.plot(np.arange(N), self.RT_path, color='blue')[0]
c1 = line1.get_color()
ax.set_xlabel('Periods')
ax.set_ylabel('Cumulative return', color=c1)
ax_ = ax.twinx()
line2 = ax_.plot(np.arange(N), self.g_path, ls='--', color='green')[0]
c2 = line2.get_color()
ax_.set_ylabel('Government expenditures', color=c2)
plt.show()
print(f"P \n {self.cp.P}")
print(f"Q \n {Q}")
print(f"Govt expenditures in {', '.join(self.states)} = {self.cp.y.flatten()}
↪")
print(f"Constant tax collections = {self.T_bar}")
print(f"Govt debt in {len(self.states)} states = {-self.b}")
print("")
print(f"Government tax collections minus debt levels in {', '.join(self.
↪states)}")
for i in range(len(self.states)):
TB = self.T_bar + self.b[i]
print(f" T+b in {self.states[i]} = {TB}")
print("")
print(f"Total government spending in {', '.join(self.states)}")
for i in range(len(self.states)):
G = self.cp.y[i, 0] + Q[i, :] @ self.b
print(f" {self.states[i]} = {G}")
print("")
print("Let's see ex-post and ex-ante returns on Arrow securities \n")
print("")
exant = 1 / self.cp.β
print(f"Ex-ante returns to purchase of Arrow securities = {exant}")
print("")
print("The Ex-post one-period gross return on the portfolio of government␣
↪assets")
print(self.R)
print("")
print("The cumulative return earned from holding 1 unit market portfolio of␣
↪government bonds")
print(self.RT_path[-1])
7.4.1 Parameters
γ = .1
λ = .1
ϕ = .1
θ = .1
ψ = .1
g_L = .5
g_M = .8
g_H = 1.2
β = .96
7.4.2 Example 1
This example is designed to produce some stylized versions of tax, debt, and deficit paths followed by the United States
during and after the Civil War and also during and after World War I.
We set the Markov chain to have three states
1−𝜆 𝜆 0
𝑃 =⎡
⎢ 0 1 − 𝜙 𝜙⎤⎥
⎣ 0 0 1⎦
P
[[0.9 0.1 0. ]
[0. 0.9 0.1]
[0. 0. 1. ]]
Q
[[0.864 0.096 0. ]
[0. 0.864 0.096]
[0. 0. 0.96 ]]
Govt expenditures in peace, war, postwar = [0.5 1.2 0.8]
Constant tax collections = 0.7548096885813149
Govt debt in 3 states = [-1. -4.07093426 -1.12975779]
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
0.17908622141460231
# The following shows the use of the wrapper class when a specific state path is given
s_path = [0, 0, 1, 1, 2]
ts_s_path = TaxSmoothingExample(g_ex1, P_ex1, b0_ex1, states_ex1, s_path=s_path)
ts_s_path.display()
P
[[0.9 0.1 0. ]
[0. 0.9 0.1]
[0. 0. 1. ]]
Q
[[0.864 0.096 0. ]
[0. 0.864 0.096]
[0. 0. 0.96 ]]
Govt expenditures in peace, war, postwar = [0.5 1.2 0.8]
Constant tax collections = 0.7548096885813149
Govt debt in 3 states = [-1. -4.07093426 -1.12975779]
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
0.9045311615620277
7.4.3 Example 2
This example captures a peace followed by a war, eventually followed by a permanent peace .
Here we set
1 0 0
𝑃 =⎡
⎢0 1−𝛾 𝛾 ⎤⎥
⎣𝜙 0 1 − 𝜙⎦
P
[[1. 0. 0. ]
[0. 0.9 0.1]
[0.1 0. 0.9]]
Q
[[0.96 0. 0. ]
[0. 0.864 0.096]
[0.096 0. 0.864]]
Govt expenditures in peace, temporary peace, war = [0.5 0.5 1.2]
Constant tax collections = 0.6053287197231834
Govt debt in 3 states = [ 2.63321799 -1. -2.51384083]
Government tax collections minus debt levels in peace, temporary peace, war
T+b in peace = -2.0278892733564
T+b in temporary peace = 1.6053287197231834
T+b in war = 3.1191695501730106
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
-9.368991732594216
7.4.4 Example 3
This example features a situation in which one of the states is a war state with no hope of peace next period, while another
state is a war state with a positive probability of peace next period.
The Markov chain is:
1−𝜆 𝜆 0 0
⎡ 0 1−𝜙 𝜙 0 ⎤
𝑃 =⎢ ⎥
⎢ 0 0 1−𝜓 𝜓 ⎥
⎣ 𝜃 0 0 1 − 𝜃⎦
with government expenditure levels for the four states being [𝑔𝐿 𝑔𝐿 𝑔𝐻 𝑔𝐻 ] where 𝑔𝐿 < 𝑔𝐻 .
We start with 𝑏0 = 1 and 𝑠0 = 1.
P
[[0.9 0.1 0. 0. ]
[0. 0.9 0.1 0. ]
[0. 0. 0.9 0.1]
[0.1 0. 0. 0.9]]
Q
[[0.864 0.096 0. 0. ]
[0. 0.864 0.096 0. ]
[0. 0. 0.864 0.096]
[0.096 0. 0. 0.864]]
Govt expenditures in peace1, peace2, war1, war2 = [0.5 0.5 1.2 1.2]
Constant tax collections = 0.6927944572748268
Govt debt in 4 states = [-1. -3.42494226 -6.86027714 -4.43533487]
Government tax collections minus debt levels in peace1, peace2, war1, war2
T+b in peace1 = 1.6927944572748268
T+b in peace2 = 4.117736720554273
T+b in war1 = 7.553071593533488
T+b in war2 = 5.128129330254041
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
0.02371440178864222
7.4.5 Example 4
with government expenditure levels for the five states being [𝑔𝐿 𝑔𝐿 𝑔𝐻 𝑔𝐻 𝑔𝐿 ] where 𝑔𝐿 < 𝑔𝐻 .
We ssume that 𝑏0 = 1 and 𝑠0 = 1.
P
[[0.9 0.1 0. 0. 0. ]
[0. 0.9 0.1 0. 0. ]
[0. 0. 0.9 0.1 0. ]
[0. 0. 0. 0.9 0.1]
[0. 0. 0. 0. 1. ]]
Q
[[0.864 0.096 0. 0. 0. ]
[0. 0.864 0.096 0. 0. ]
[0. 0. 0.864 0.096 0. ]
[0. 0. 0. 0.864 0.096]
[0. 0. 0. 0. 0.96 ]]
Govt expenditures in peace1, peace2, war1, war2, permanent peace = [0.5 0.5 1.2 1.
↪2 0.5]
Government tax collections minus debt levels in peace1, peace2, war1, war2,␣
↪permanent peace
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
-11.132109773063616
7.4.6 Example 5
The example captures a case when the system follows a deterministic path from peace to war, and back to peace again.
Since there is no randomness, the outcomes in complete markets setting should be the same as in incomplete markets
setting.
The Markov chain is:
0 1 0 0 0 0 0
⎡0 0 1 0 0 0 0⎤
⎢ ⎥
⎢0 0 0 1 0 0 0⎥
𝑃 = ⎢0 0 0 0 1 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 1 0⎥
⎢0 0 0 0 0 0 1⎥
⎣0 0 0 0 0 0 1⎦
with government expenditure levels for the seven states being [𝑔𝐿 𝑔𝐿 𝑔𝐻 𝑔𝐻 𝑔𝐻 𝑔𝐻 𝑔𝐿 ] where 𝑔𝐿 < 𝑔𝐻 .
Assume 𝑏0 = 1 and 𝑠0 = 1.
ts_ex5.display()
P
[[0 1 0 0 0 0 0]
[0 0 1 0 0 0 0]
[0 0 0 1 0 0 0]
[0 0 0 0 1 0 0]
[0 0 0 0 0 1 0]
[0 0 0 0 0 0 1]
[0 0 0 0 0 0 1]]
Q
[[0. 0.96 0. 0. 0. 0. 0. ]
[0. 0. 0.96 0. 0. 0. 0. ]
[0. 0. 0. 0.96 0. 0. 0. ]
[0. 0. 0. 0. 0.96 0. 0. ]
[0. 0. 0. 0. 0. 0.96 0. ]
[0. 0. 0. 0. 0. 0. 0.96]
[0. 0. 0. 0. 0. 0. 0.96]]
Govt expenditures in peace1, peace2, war1, war2, war3, permanent peace = [0.5 0.5␣
↪1.2 1.2 1.2 1.2 0.5]
Government tax collections minus debt levels in peace1, peace2, war1, war2, war3,␣
↪permanent peace
Total government spending in peace1, peace2, war1, war2, war3, permanent peace
peace1 = 1.5571895472128
peace2 = 1.6584286588928003
war1 = 1.7638860668928
war2 = 1.1445708668928003
war3 = 0.49945086689280027
permanent peace = -0.17254913310719933
The cumulative return earned from holding 1 unit market portfolio of government␣
↪bonds
1.2775343959060068
To construct a tax-smoothing version of the complete markets consumption-smoothing model with a continuous state
space that we presented in the lecture consumption smoothing with complete and incomplete markets, we simply relabel
variables.
Thus, a government faces a sequence of budget constraints
where 𝑇𝑡 is tax revenues, 𝑏𝑡 are receipts at 𝑡 from contingent claims that the government had purchased at time 𝑡 − 1, and
which states that the present value of government purchases equals the value of government assets at 𝑡 plus the present
value of tax receipts.
With these relabelings, examples presented in consumption smoothing with complete and incomplete markets can be inter-
preted as tax-smoothing models.
Returns: In the continuous state version of our incomplete markets model, the ex post one-period gross rate of return
on the government portfolio equals
𝑏(𝑥𝑡+1 )
𝑅(𝑥𝑡+1 |𝑥𝑡 ) =
𝛽𝐸𝑏(𝑥𝑡+1 )|𝑥𝑡
Related Lectures
Throughout this lecture, we have taken one-period interest rates and Arrow security prices as exogenous objects deter-
mined outside the model and specified them in ways designed to align our models closely with the consumption smoothing
model of Barro [Barro, 1979].
Other lectures make these objects endogenous and describe how a government optimally manipulates prices of govern-
ment debt, albeit indirectly via effects distorting taxes have on equilibrium prices and allocations.
In optimal taxation in an LQ economy and recursive optimal taxation, we study complete-markets models in which the
government recognizes that it can manipulate Arrow securities prices.
Linear-quadratic versions of the Lucas-Stokey tax-smoothing model are described in Optimal Taxation in an LQ Economy.
That lecture is a warm-up for the non-linear-quadratic model of tax smoothing described in Optimal Taxation with State-
Contingent Debt.
In both Optimal Taxation in an LQ Economy and Optimal Taxation with State-Contingent Debt, the government recognizes
that its decisions affect prices.
In optimal taxation with incomplete markets, we study an incomplete-markets model in which the government also
manipulates prices of government debt.
EIGHT
In addition to what’s in Anaconda, this lecture will need the following libraries:
8.1 Overview
This lecture describes Markov jump linear quadratic dynamic programming, an extension of the method described
in the first LQ control lecture.
Markov jump linear quadratic dynamic programming is described and analyzed in [Do Val et al., 1999] and the references
cited there.
The method has been applied to problems in macroeconomics and monetary economics by [Svensson et al., 2008] and
[Svensson and Williams, 2009].
The periodic models of seasonality described in chapter 14 of [Hansen and Sargent, 2013] are a special case of Markov
jump linear quadratic problems.
Markov jump linear quadratic dynamic programming combines advantages of
• the computational simplicity of linear quadratic dynamic programming, with
• the ability of finite state Markov chains to represent interesting patterns of random variation.
The idea is to replace the constant matrices that define a linear quadratic dynamic programming problem with 𝑁 sets
of matrices that are fixed functions of the state of an 𝑁 state Markov chain.
The state of the Markov chain together with the continuous 𝑛 × 1 state vector 𝑥𝑡 form the state of the system.
For the class of infinite horizon problems being studied in this lecture, we obtain 𝑁 interrelated matrix Riccati equations
that determine 𝑁 optimal value functions and 𝑁 linear decision rules.
One of these value functions and one of these decision rules apply in each of the 𝑁 Markov states.
That is, when the Markov state is in state 𝑗, the value function and the decision rule for state 𝑗 prevails.
143
Advanced Quantitative Economics with Python
𝑢𝑡 = −𝐹 𝑥𝑡
− (𝑥′𝑡 𝑃 𝑥𝑡 + 𝜌)
𝜌 = 𝛽 (𝜌 + trace(𝑃 𝐶𝐶 ′ ))
With the preceding formulas in mind, we are ready to approach Markov Jump linear quadratic dynamic programming.
The key idea is to make the matrices 𝐴, 𝐵, 𝐶, 𝑅, 𝑄, 𝑊 fixed functions of a finite state 𝑠 that is governed by an 𝑁 state
Markov chain.
This makes decision rules depend on the Markov state, and so fluctuate through time in limited ways.
In particular, we use the following extension of a discrete-time linear quadratic dynamic programming problem.
We let 𝑠𝑡 ∈ [1, 2, … , 𝑁 ] be a time 𝑡 realization of an 𝑁 -state Markov chain with transition matrix Π having typical
element Π𝑖𝑗 .
We’ll switch between labeling today’s state as 𝑠𝑡 and 𝑖 and between labeling tomorrow’s state as 𝑠𝑡+1 or 𝑗.
The decision-maker solves the minimization problem:
∞
min
∞
𝐸 ∑ 𝛽 𝑡 𝑟(𝑥𝑡 , 𝑠𝑡 , 𝑢𝑡 )
{𝑢𝑡 }𝑡=0
𝑡=0
with
subject to linear laws of motion with matrices (𝐴, 𝐵, 𝐶) each possibly dependent on the Markov-state-𝑠𝑡 :
𝑢𝑡 = −𝐹𝑠𝑡 𝑥𝑡
or equivalently
−𝑥′𝑡 𝑃𝑖 𝑥𝑡 − 𝜌𝑖
The optimal value functions −𝑥′ 𝑃𝑖 𝑥 − 𝜌𝑖 for 𝑖 = 1, … , 𝑛 satisfy the 𝑁 interrelated Bellman equations
−𝑥′ 𝑃𝑖 𝑥 − 𝜌𝑖 = max −𝑥′ 𝑅𝑖 𝑥 +𝑢′ 𝑄𝑖 𝑢 + 2𝑢′ 𝑊𝑖 𝑥 + 𝛽 ∑ Π𝑖𝑗 𝐸((𝐴𝑖 𝑥 + 𝐵𝑖 𝑢 + 𝐶𝑖 𝑤)′ 𝑃𝑗 (𝐴𝑖 𝑥 + 𝐵𝑖 𝑢 + 𝐶𝑖 𝑤)𝑥 + 𝜌𝑗 )
𝑢
𝑗
The matrices 𝑃𝑠𝑡 = 𝑃𝑖 and the scalars 𝜌𝑠𝑡 = 𝜌𝑖 , 𝑖 = 1, …, n satisfy the following stacked system of algebraic matrix
Riccati equations:
8.4 Applications
import numpy as np
import quantecon as qe
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
8.5 Example 1
This example is a version of a classic problem of optimally adjusting a variable 𝑘𝑡 to a target level in the face of costly
adjustment.
This provides a model of gradual adjustment.
Given 𝑘0 , the objective function is
∞
max
∞
𝐸0 ∑ 𝛽 𝑡 𝑟 (𝑠𝑡 , 𝑘𝑡 )
{𝑘𝑡 }𝑡=1 𝑡=0
𝐸0 is a mathematical expectation conditioned on time 0 information 𝑥0 , 𝑠0 and the transition law for continuous state
variable 𝑘𝑡 is
𝑘𝑡+1 − 𝑘𝑡 = 𝑢𝑡
We can think of 𝑘𝑡 as the decision-maker’s capital and 𝑢𝑡 as costs of adjusting the level of capital.
We assume that 𝑓1 (𝑠𝑡 ) > 0, 𝑓2 (𝑠𝑡 ) > 0, and 𝑑 (𝑠𝑡 ) > 0.
Denote the state transition matrix for Markov state 𝑠𝑡 ∈ {1, 2} as Π:
Pr (𝑠𝑡+1 = 𝑗 ∣ 𝑠𝑡 = 𝑖) = Π𝑖𝑗
𝑘
Let 𝑥𝑡 = [ 𝑡 ]
1
We can represent the one-period payoff function 𝑟 (𝑠𝑡 , 𝑘𝑡 ) as
Rs = np.zeros((m, n, n))
Qs = np.zeros((m, k, k))
for i in range(m):
Rs[i, 0, 0] = f2_vals[i]
Rs[i, 1, 0] = - f1_vals[i] / 2
Rs[i, 0, 1] = - f1_vals[i] / 2
Qs[i, 0, 0] = d_vals[i]
The continuous part of the state 𝑥𝑡 consists of two variables, namely, 𝑘𝑡 and a constant term.
We start with a Markov transition matrix that makes the Markov state be strictly periodic:
0 1
Π1 = [ ],
1 0
We set 𝑓1,𝑠𝑡 and 𝑓2,𝑠𝑡 to be independent of the Markov state 𝑠𝑡
𝑓1,1 = 𝑓1,2 = 1,
𝑓2,1 = 𝑓2,2 = 1
In contrast to 𝑓1,𝑠𝑡 and 𝑓2,𝑠𝑡 , we make the adjustment cost 𝑑𝑠𝑡 vary across Markov states 𝑠𝑡 .
We set the adjustment cost to be lower in Markov state 2
𝑑1 = 1, 𝑑2 = 0.5
The following code forms a Markov switching LQ problem and computes the optimal value functions and optimal decision
rules for each Markov state
# Construct matrices
Qs, Rs, Ns, As, Bs, Cs, k_star = construct_arrays1(d_vals=[1., 0.5])
Let’s look at the value function matrices and the decision rules for each Markov state
# P(s)
ex1_a.Ps
[[ 1.37424214, -0.68712107],
[-0.68712107, -4.65643947]]])
array([0., 0.])
# F(s)
ex1_a.Fs
[[ 0.74848427, -0.37424214]]])
Now we’ll plot the decision rules and see if they make sense
fig, ax = plt.subplots()
ax.plot(k_grid, k_grid + u1_star, label=r"$\overline{s}_1$ (high)")
ax.plot(k_grid, k_grid + u2_star, label=r"$\overline{s}_2$ (low)")
# The optimal k*
ax.scatter([0.5, 0.5], [0.5, 0.5], marker="*")
ax.plot([k_star[0], k_star[0]], [0., 1.0], '--')
(continues on next page)
# 45 degree line
ax.plot([0., 1.], [0., 1.], '--', label="45 degree line")
ax.set_xlabel("$k_t$")
ax.set_ylabel("$k_{t+1}$")
ax.legend()
plt.show()
The above graph plots 𝑘𝑡+1 = 𝑘𝑡 + 𝑢𝑡 = 𝑘𝑡 − 𝐹 𝑥𝑡 as an affine (i.e., linear in 𝑘𝑡 plus a constant) function of 𝑘𝑡 for both
Markov states 𝑠𝑡 .
It also plots the 45 degree line.
Notice that the two 𝑠𝑡 -dependent closed loop functions that determine 𝑘𝑡+1 as functions of 𝑘𝑡 share the same rest point
(also called a fixed point) at 𝑘𝑡 = 0.5.
Evidently, the optimal decision rule in Markov state 2, in which the adjustment cost is lower, makes 𝑘𝑡+1 a flatter function
of 𝑘𝑡 in Markov state 2.
This happens because when 𝑘𝑡 is not at its fixed point, |𝑢𝑡,2 | > |𝑢𝑡,2 |, so that the decision-maker adjusts toward the fixed
point faster when the Markov state 𝑠𝑡 takes a value that makes it cheaper.
Now we’ll depart from the preceding transition matrix that made the Markov state be strictly periodic.
We’ll begin with symmetric transition matrices of the form
1−𝜆 𝜆
Π2 = [ ].
𝜆 1−𝜆
λ = 0.8 # high λ
Π2 = np.array([[1-λ, λ],
[λ, 1-λ]])
[[ 0.74434525, -0.37217263]]])
λ = 0.2 # low λ
Π2 = np.array([[1-λ, λ],
[λ, 1-λ]])
[[ 0.72818728, -0.36409364]]])
for i, λ in enumerate(λ_vals):
Π2 = np.array([[1-λ, λ],
[λ, 1-λ]])
ax.set_xlabel(r"$\lambda$")
ax.set_ylabel("$F_{s_t}$")
ax.set_title(f"Coefficient on {state_var}")
ax.legend()
plt.show()
Notice how the decision rules’ constants and slopes behave as functions of 𝜆.
Evidently, as the Markov chain becomes more nearly periodic (i.e., as 𝜆 → 1), the dynamic program adjusts capital faster
in the low adjustment cost Markov state to take advantage of what is only temporarily a more favorable time to invest.
Now let’s study situations in which the Markov transition matrix Π is asymmetric
1−𝜆 𝜆
Π3 = [ ].
𝛿 1−𝛿
λ, δ = 0.8, 0.2
Π3 = np.array([[1-λ, λ],
[δ, 1-δ]])
[[ 0.72749075, -0.36374537]]])
for i, λ in enumerate(λ_vals):
λ_grid[i, :] = λ
δ_grid[i, :] = δ_vals
for j, δ in enumerate(δ_vals):
Π3 = np.array([[1-λ, λ],
[δ, 1-δ]])
The following code defines a wrapper function that computes optimal decision rules for cases with different Markov
transition matrices
# Symmetric Π
# Notice that pure periodic transition is a special case
# when λ=1
print("symmetric Π case:\n")
λ_vals = np.linspace(0., 1., 10)
F1 = np.empty((λ_vals.size, len(state_vec)))
F2 = np.empty((λ_vals.size, len(state_vec)))
for i, λ in enumerate(λ_vals):
Π2 = np.array([[1-λ, λ],
[λ, 1-λ]])
ax.set_xlabel(r"$\lambda$")
ax.set_ylabel(r"$F(\overline{s}_t)$")
ax.set_title(f"coefficient on {state_var}")
ax.legend()
plt.show()
ax.set_xlabel(r"$\lambda$")
ax.set_ylabel("$k$")
ax.set_title("Optimal k levels and k targets")
ax.text(0.5, min(k_star)+(max(k_star)-min(k_star))/20, r"$\lambda=0.5$")
(continues on next page)
# Asymmetric Π
print("asymmetric Π case:\n")
δ_vals = np.linspace(0., 1., 10)
for i, λ in enumerate(λ_vals):
λ_grid[i, :] = λ
δ_grid[i, :] = δ_vals
for j, δ in enumerate(δ_vals):
Π3 = np.array([[1-λ, λ],
[δ, 1-δ]])
To illustrate the code with another example, we shall set 𝑓2,𝑠𝑡 and 𝑑𝑠𝑡 as constant functions and
Thus, the sole role of the Markov jump state 𝑠𝑡 is to identify times in which capital is very productive and other times in
which it is less productive.
The example below reveals much about the structure of the optimum problem and optimal policies.
Only 𝑓1,𝑠𝑡 varies with 𝑠𝑡 .
𝑓1,𝑠𝑡
So there are different 𝑠𝑡 -dependent optimal static 𝑘 level in different states 𝑘𝑠∗𝑡 = 2𝑓2,𝑠𝑡 , values of 𝑘 that maximize
one-period payoff functions in each state.
We denote a target 𝑘 level as 𝑘𝑠𝑡𝑎𝑟𝑔𝑒𝑡
𝑡
, the fixed point of the optimal policies in each state, given the value of 𝜆.
We call 𝑘𝑠𝑡𝑎𝑟𝑔𝑒𝑡
𝑡
a “target” because in each Markov state 𝑠𝑡 , optimal policies are contraction mappings and will push 𝑘𝑡
towards a fixed point 𝑘𝑠𝑡𝑎𝑟𝑔𝑒𝑡
𝑡
.
When 𝜆 → 0, each Markov state becomes close to absorbing state and consequently 𝑘𝑠𝑡𝑎𝑟𝑔𝑒𝑡
𝑡
→ 𝑘𝑠∗𝑡 .
But when 𝜆 → 1, the Markov transition matrix becomes more nearly periodic, so the optimum decision rules target more
at the optimal 𝑘 level in the other state in order to enjoy higher expected payoff in the next period.
The switch happens at 𝜆 = 0.5 when both states are equally likely to be reached.
Below we plot an additional figure that shows optimal 𝑘 levels in the two states Markov jump state and also how the
targeted 𝑘 levels change as 𝜆 changes.
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
8.6 Example 2
We now add to the example 1 setup another state variable 𝑤𝑡 that follows the evolution law
We think of 𝑤𝑡 as a rental rate or tax rate that the decision maker pays each period for 𝑘𝑡 .
To capture this idea, we add to the decision-maker’s one-period payoff function the product of 𝑤𝑡 and 𝑘𝑡
𝑘𝑡
We now let the continuous part of the state at time 𝑡 be 𝑥𝑡 = ⎡ ⎤
⎢ 1 ⎥ and continue to set the control 𝑢𝑡 = 𝑘𝑡+1 − 𝑘𝑡 .
⎣𝑤𝑡 ⎦
We can write the one-period payoff function 𝑟 (𝑠𝑡 , 𝑘𝑡 , 𝑤𝑡 ) as
2
𝑟 (𝑠𝑡 , 𝑘𝑡 , 𝑤𝑡 ) = 𝑓1 (𝑠𝑡 ) 𝑘𝑡 − 𝑓2 (𝑠𝑡 ) 𝑘𝑡2 − 𝑑 (𝑠𝑡 ) (𝑘𝑡+1 − 𝑘𝑡 ) − 𝑤𝑡 𝑘𝑡
⎛ ⎞
⎜
⎜ 𝑓2 (𝑠𝑡 ) − 𝑓1 (𝑠 2
𝑡) 1
2
⎟
⎟
⎜ ′ ⎡ 𝑓 (𝑠 ) ⎤ 2⎟
= −⎜
⎜𝑥 𝑡 ⎢− 1
2
𝑡
0 0 ⎥𝑥𝑡 + 𝑑
⏟ (𝑠 𝑡 ) 𝑢𝑡 ⎟
⎟ ,
⎜
⎜ ⏟⏟ 1
0 0 ≡𝑄(𝑠𝑡 ) ⎟
⎟
⎣ ⏟ 2⏟⏟⏟⏟⏟⏟⏟⏟ ⎦
⎝ ≡𝑅(𝑠𝑡 ) ⎠
m = len(f1_vals)
n, k, j = 3, 1, 1
Rs = np.zeros((m, n, n))
Qs = np.zeros((m, k, k))
As = np.zeros((m, n, n))
Bs = np.zeros((m, n, k))
Cs = np.zeros((m, n, j))
for i in range(m):
Rs[i, 0, 0] = f2_vals[i]
Rs[i, 1, 0] = - f1_vals[i] / 2
(continues on next page)
Qs[i, 0, 0] = d_vals[i]
As[i, 0, 0] = 1
As[i, 1, 1] = 1
As[i, 2, 1] = α0_vals[i]
As[i, 2, 2] = ρ_vals[i]
Ns = None
k_star = None
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
symmetric Π case:
asymmetric Π case:
The following lectures describe how Markov jump linear quadratic dynamic programming can be used to extend the
[Barro, 1979] model of optimal tax-smoothing and government debt in several interesting directions
1. How to Pay for a War: Part 1
2. How to Pay for a War: Part 2
3. How to Pay for a War: Part 3
NINE
9.1 Overview
This lecture uses the method of Markov jump linear quadratic dynamic programming that is described in lecture
Markov Jump LQ dynamic programming to extend the [Barro, 1979] model of optimal tax-smoothing and government
debt in a particular direction.
This lecture has two sequels that offer further extensions of the Barro model
1. How to Pay for a War: Part 2
2. How to Pay for a War: Part 3
The extensions are modified versions of his 1979 model suggested by [Barro, 1999] and [Barro and McCleary, 2003]).
[Barro, 1979] m is about a government that borrows and lends in order to minimize an intertemporal measure of distortions
caused by taxes.
Technical tractability induced [Barro, 1979] to assume that
• the government trades only one-period risk-free debt, and
• the one-period risk-free interest rate is constant
By using Markov jump linear quadratic dynamic programming we can allow interest rates to move over time in empirically
interesting ways.
Also, by expanding the dimension of the state, we can add a maturity composition decision to the government’s problem.
By doing these two things we extend [Barro, 1979] along lines he suggested in [Barro, 1999] and [Barro and McCleary,
2003]).
[Barro, 1979] assumed
• that a government faces an exogenous sequence of expenditures that it must finance by a tax collection sequence
whose expected present value equals the initial debt it owes plus the expected present value of those expenditures.
∞
• that the government wants to minimize a measure of tax distortions that is proportional to 𝐸0 ∑𝑡=0 𝛽 𝑡 𝑇𝑡2 , where
𝑇𝑡 are total tax collections and 𝐸0 is a mathematical expectation conditioned on time 0 information.
• that the government trades only one asset, a risk-free one-period bond.
• that the gross interest rate on the one-period bond is constant and equal to 𝛽 −1 , the reciprocal of the factor 𝛽 at
which the government discounts future tax distortions.
Barro’s model can be mapped into a discounted linear quadratic dynamic programming problem.
Partly inspired by [Barro, 1999] and [Barro and McCleary, 2003], our generalizations of [Barro, 1979], assume
• that the government borrows or saves in the form of risk-free bonds of maturities 1, 2, … , 𝐻.
197
Advanced Quantitative Economics with Python
• that interest rates on those bonds are time-varying and in particular, governed by a jointly stationary stochastic
process.
Our generalizations are designed to fit within a generalization of an ordinary linear quadratic dynamic programming
problem in which matrices that define the quadratic objective function and the state transition function are time-varying
and stochastic.
This generalization, known as a Markov jump linear quadratic dynamic program, combines
• the computational simplicity of linear quadratic dynamic programming, and
• the ability of finite state Markov chains to represent interesting patterns of random variation.
We want the stochastic time variation in the matrices defining the dynamic programming problem to represent variation
over time in
• interest rates
• default rates
• roll over risks
As described in Markov Jump LQ dynamic programming, the idea underlying Markov jump linear quadratic dynamic
programming is to replace the constant matrices defining a linear quadratic dynamic programming problem with
matrices that are fixed functions of an 𝑁 state Markov chain.
For infinite horizon problems, this leads to 𝑁 interrelated matrix Riccati equations that pin down 𝑁 value functions and
𝑁 linear decision rules, applying to the 𝑁 Markov states.
import quantecon as qe
import numpy as np
import matplotlib.pyplot as plt
We begin by solving a version of [Barro, 1979] by mapping it into the original LQ framework.
As mentioned in this lecture, the Barro model is mathematically isomorphic with the LQ permanent income model.
Let
• 𝑇𝑡 denote tax collections
• 𝛽 be a discount factor
• 𝑏𝑡,𝑡+1 be time 𝑡 + 1 goods that at 𝑡 the government promises to deliver to time 𝑡 buyers of one-period government
bonds
• 𝐺𝑡 be government purchases
• 𝑝𝑡,𝑡+1 the number of time 𝑡 goods received per time 𝑡 + 1 goods promised to one-period bond purchasers.
Evidently, 𝑝𝑡,𝑡+1 is inversely related to appropriate corresponding gross interest rates on government debt.
In the spirit of [Barro, 1979], the stochastic process of government expenditures is exogenous.
The government’s problem is to choose a plan for taxation and borrowing {𝑏𝑡+1 , 𝑇𝑡 }∞
𝑡=0 to minimize
∞
𝐸0 ∑ 𝛽 𝑡 𝑇𝑡2
𝑡=0
𝐺𝑡 = 𝑈𝑔 𝑧𝑡
𝑧𝑡+1 = 𝐴22 𝑧𝑡 + 𝐶2 𝑤𝑡+1
where 𝑤𝑡+1 ∼ 𝑁 (0, 𝐼)
The variables 𝑇𝑡 , 𝑏𝑡,𝑡+1 are control variables chosen at 𝑡, while 𝑏𝑡−1,𝑡 is an endogenous state variable inherited from the
past at time 𝑡 and 𝑝𝑡,𝑡+1 is an exogenous state variable at time 𝑡.
To begin, we assume that 𝑝𝑡,𝑡+1 is constant (and equal to 𝛽)
• later we will extend the model to allow 𝑝𝑡,𝑡+1 to vary over time
𝑏𝑡−1,𝑡
To map into the LQ framework, we use 𝑥𝑡 = [ ] as the state vector, and 𝑢𝑡 = 𝑏𝑡,𝑡+1 as the control variable.
𝑧𝑡
Therefore, the (𝐴, 𝐵, 𝐶) matrices are defined by the state-transition law:
0 0 1 0
𝑥𝑡+1 = [ ] 𝑥 + [ ] 𝑢𝑡 + [ ] 𝑤𝑡+1
0 𝐴22 𝑡 0 𝐶2
To find the appropriate (𝑅, 𝑄, 𝑊 ) matrices, we note that 𝐺𝑡 and 𝑏𝑡−1,𝑡 can be written as appropriately defined functions
of the current state:
𝐺𝑡 = 𝑆𝐺 𝑥𝑡 , 𝑏𝑡−1,𝑡 = 𝑆1 𝑥𝑡
If we define 𝑀𝑡 = −𝑝𝑡,𝑡+1 , and let 𝑆 = 𝑆𝐺 + 𝑆1 , then we can write taxation as a function of the states and control using
the government’s budget constraint:
𝑇𝑡 = 𝑆𝑥𝑡 + 𝑀𝑡 𝑢𝑡
1
To do this, we set 𝑧𝑡 = [ ], and consequently:
𝐺𝑡
1 0 0
𝐴22 = [ ̄ ] , 𝐶2 = [ ]
𝐺 𝜌 𝜎
# Model parameters
β, Gbar, ρ, σ = 0.95, 5, 0.8, 1
C2 = np.array([[0],
[σ]])
Ug = np.array([[0, 1]])
# LQ framework matrices
A_t = np.zeros((1, 3))
A_b = np.hstack((np.zeros((2, 1)), A22))
A = np.vstack((A_t, A_b))
B = np.zeros((3, 1))
B[0, 0] = 1
M = np.array([[-β]])
R = S.T @ S
Q = M.T @ M
W = M.T @ S
We can see the isomorphism by noting that consumption is a martingale in the permanent income model and that taxation
is a martingale in Barro’s model.
We can check this using the 𝐹 matrix of the LQ model.
Because 𝑢𝑡 = −𝐹 𝑥𝑡 , we have
𝑇𝑡 = 𝑆𝑥𝑡 + 𝑀 𝑢𝑡 = (𝑆 − 𝑀 𝐹 )𝑥𝑡
and
(𝑆 − 𝑀 𝐹 )(𝐴 − 𝐵𝐹 ) = (𝑆 − 𝑀 𝐹 ),
S - M @ F, (S - M @ F) @ (A - B @ F)
This explains the fanning out of the conditional empirical distribution of taxation across time, computed by simulating
the Barro model many times and averaging over simulated paths:
T = 500
for i in range(250):
x, u, w = LQBarro.compute_sequence(x0, ts_length=T)
plt.plot(list(range(T+1)), ((S - M @ F) @ x)[0, :])
plt.xlabel('Time')
plt.ylabel('Taxation')
plt.show()
We can see a similar, but a smoother pattern, if we plot government debt over time.
T = 500
for i in range(250):
x, u, w = LQBarro.compute_sequence(x0, ts_length=T)
plt.plot(list(range(T+1)), x[0, :])
plt.xlabel('Time')
plt.ylabel('Govt Debt')
plt.show()
To implement the extension to the Barro model in which 𝑝𝑡,𝑡+1 varies over time, we must allow the M matrix to be
time-varying.
Our 𝑄 and 𝑊 matrices must also vary over time.
We can solve such a model using the LQMarkov class that solves Markov jump linear quandratic control problems as
described above.
The code for the class can be viewed here.
The class takes lists of matrices that corresponds to 𝑁 Markov states.
The value and policy functions are then found by iterating on a coupled system of matrix Riccati difference equations.
Optimal 𝑃𝑠 , 𝐹𝑠 , 𝑑𝑠 are stored as attributes.
The class also contains a method that simulates a model.
9.4. Python Class to Solve Markov Jump Linear Quadratic Control Problems 203
Advanced Quantitative Economics with Python
We can use the above class to implement a version of the Barro model with a time-varying interest rate.
A simple way to extend the model is to allow the interest rate to take two possible values.
We set:
1
𝑝𝑡,𝑡+1 = 𝛽 + 0.02 = 0.97
2
𝑝𝑡,𝑡+1 = 𝛽 − 0.017 = 0.933
Thus, the first Markov state has a low interest rate and the second Markov state has a high interest rate.
We must also specify a transition matrix for the Markov state.
We use:
0.8 0.2
Π=[ ]
0.2 0.8
Here, each Markov state is persistent, and there is are equal chances of moving from one state to the other.
The choice of parameters means that the unconditional expectation of 𝑝𝑡,𝑡+1 is 0.9515, higher than 𝛽(= 0.95).
If we were to set 𝑝𝑡,𝑡+1 = 0.9515 in the version of the model with a constant interest rate, government debt would
explode.
As = [A, A]
Bs = [B, B]
Cs = [C, C]
Rs = [R, R]
M1 = np.array([[-β - 0.02]])
M2 = np.array([[-β + 0.017]])
Q1 = M1.T @ M1
Q2 = M2.T @ M2
Qs = [Q1, Q2]
W1 = M1.T @ S
W2 = M2.T @ S
Ws = [W1, W2]
lqm.Fs[0]
lqm.Fs[1]
Simulating a large number of such economies over time reveals interesting dynamics.
Debt tends to stay low and stable but recurrently surges.
T = 2000
x0 = np.array([[1000, 1, 25]])
for i in range(250):
x, u, w, s = lqm.compute_sequence(x0, ts_length=T)
plt.plot(list(range(T+1)), x[0, :])
plt.xlabel('Time')
plt.ylabel('Govt Debt')
plt.show()
TEN
10.1 Overview
This lecture presents another application of Markov jump linear quadratic dynamic programming and constitutes a sequel
to an earlier lecture.
We use a method introduced in lecture Markov Jump LQ dynamic programming toimplement suggestions by [Barro, 1999]
and [Barro and McCleary, 2003]) for extending his classic 1979 model of tax smoothing.
[Barro, 1979] model is about a government that borrows and lends in order to help it minimize an intertemporal measure
of distortions caused by taxes.
Technically, [Barro, 1979] model looks a lot like a consumption-smoothing model.
Our generalizations of [Barro, 1979] will also look like souped-up consumption-smoothing models.
Wanting tractability induced [Barro, 1979] to assume that
• the government trades only one-period risk-free debt, and
• the one-period risk-free interest rate is constant
In our earlier lecture, we relaxed the second of these assumptions but not the first.
In particular, we used Markov jump linear quadratic dynamic programming to allow the exogenous interest rate to vary
over time.
In this lecture, we add a maturity composition decision to the government’s problem by expanding the dimension of the
state.
We assume
• that the government borrows or saves in the form of risk-free bonds of maturities 1, 2, … , 𝐻.
• that interest rates on those bonds are time-varying and in particular are governed by a jointly stationary stochastic
process.
In addition to what’s in Anaconda, this lecture deploys the quantecon library:
import quantecon as qe
import numpy as np
import matplotlib.pyplot as plt
207
Advanced Quantitative Economics with Python
Let
• 𝑇𝑡 denote tax collections
• 𝛽 be a discount factor
• 𝑏𝑡,𝑡+1 be time 𝑡 + 1 goods that the government promises to pay at 𝑡
• 𝑏𝑡,𝑡+2 betime 𝑡 + 2 goods that the government promises to pay at time 𝑡
• 𝐺𝑡 be government purchases
• 𝑝𝑡,𝑡+1 be the number of time 𝑡 goods received per time 𝑡 + 1 goods promised
• 𝑝𝑡,𝑡+2 be the number of time 𝑡 goods received per time 𝑡 + 2 goods promised.
Evidently, 𝑝𝑡,𝑡+1 , 𝑝𝑡,𝑡+2 are inversely related to appropriate corresponding gross interest rates on government debt.
In the spirit of [Barro, 1979], government expenditures are governed by an exogenous stochastic process.
Given initial conditions 𝑏−2,0 , 𝑏−1,0 , 𝑧0 , 𝑖0 , where 𝑖0 is the initial Markov state, the government chooses a contingency
plan for {𝑏𝑡,𝑡+1 , 𝑏𝑡,𝑡+2 , 𝑇𝑡 }∞
𝑡=0 to maximize.
∞
−𝐸0 ∑ 𝛽 𝑡 [𝑇𝑡2 + 𝑐1 (𝑏𝑡,𝑡+1 − 𝑏𝑡,𝑡+2 )2 ]
𝑡=0
Here
• 𝑤𝑡+1 ∼ 𝑁 (0, 𝐼) and Π𝑖𝑗 is the probability that the Markov state moves from state 𝑖 to state 𝑗 in one period
• 𝑇𝑡 , 𝑏𝑡,𝑡+1 , 𝑏𝑡,𝑡+2 are control variables chosen at time 𝑡
• variables 𝑏𝑡−1,𝑡 , 𝑏𝑡−2,𝑡 are endogenous state variables inherited from the past at time 𝑡
• 𝑝𝑡,𝑡+1 , 𝑝𝑡,𝑡+2 are exogenous state variables at time 𝑡
The parameter 𝑐1 imposes a penalty on the government’s issuing different quantities of one and two-period debt.
This penalty deters the government from taking large “long-short” positions in debt of different maturities.
An example below will show the penalty in action.
As well as extending the model to allow for a maturity decision for government debt, we can also in principle allow the
matrices 𝑈𝑔,𝑠𝑡 , 𝐴22,𝑠𝑡 , 𝐶2,𝑠𝑡 to depend on the Markov state 𝑠𝑡 .
Below, we will often adopt the convention that for matrices appearing in a linear state space, 𝐴𝑡 ≡ 𝐴𝑠𝑡 , 𝐶𝑡 ≡ 𝐶𝑠𝑡 and
so on, so that dependence on 𝑡 is always intermediated through the Markov state 𝑠𝑡 .
First, define
𝑏̂𝑡
𝑏̄𝑡 = [ ]
𝑏𝑡−1,𝑡+1
𝑏̄
𝑥𝑡 = [ 𝑡 ]
𝑧𝑡
𝑏̂ 0 1 𝑏̂𝑡 1 0 𝑏𝑡,𝑡+1
[ 𝑡+1 ] = [ ][ ]+[ ][ ]
𝑏𝑡,𝑡+2 0 0 𝑏𝑡−1,𝑡+1 0 1 𝑏𝑡,𝑡+2
or
𝐺𝑡 = 𝑆𝐺,𝑡 𝑥𝑡 , 𝑏̂𝑡 = 𝑆1 𝑥𝑡
and
𝑀𝑡 = [−𝑝𝑡,𝑡+1 −𝑝𝑡,𝑡+2 ]
where 𝑝𝑡,𝑡+1 is the discount on one period loans in the discrete Markov state at time 𝑡 and 𝑝𝑡,𝑡+2 is the discount on
two-period loans in the discrete Markov state.
Define
𝑆𝑡 = 𝑆𝐺,𝑡 + 𝑆1
𝑇𝑡 = 𝑀𝑡 𝑢𝑡 + 𝑆𝑡 𝑥𝑡
It follows that
or
where
Because the payoff function also includes the penalty parameter on issuing debt of different maturities, we have:
1 −1
where 𝑄𝑐 = [ ].
−1 1
Therefore, the appropriate 𝑄 matrix in the Markov jump LQ problem is:
𝑄𝑐𝑡 = 𝑄𝑡 + 𝑐1 𝑄𝑐
where
𝐴11 0 𝐵1 0
𝐴𝑡 = [ ], 𝐵=[ ], 𝐶𝑡 = [ ]
0 𝐴22,𝑡 0 𝐶2,𝑡
Thus, in this problem all the matrices apart from 𝐵 may depend on the Markov state at time 𝑡.
As shown in the previous lecture, when provided with appropriate 𝐴, 𝐵, 𝐶, 𝑅, 𝑄, 𝑊 matrices for each Markov state the
LQMarkov class can solve Markov jump LQ problems.
The function below maps the primitive matrices and parameters from the above two-period model into the matrices that
the LQMarkov class requires:
"""
Function which takes A22, C2, Ug, p_{t, t+1}, p_{t, t+2} and penalty
parameter c1, and returns the required matrices for the LQMarkov
model: A, B, C, R, Q, W.
This version uses the condensed version of the endogenous state.
"""
B1 = np.eye(2)
# Create M matrix
M = np.hstack((-p1, -p2))
# Create A, B, C matrices
A_T = np.hstack((A11, np.zeros((2, nz))))
A_B = np.hstack((np.zeros((nz, 2)), A22))
A = np.vstack((A_T, A_B))
# Create R, Q, W matrices
R = S.T @ S
Q = M.T @ M + c1 * Qc
W = M.T @ S
return A, B, C, R, Q, W
With the above function, we can proceed to solve the model in two steps:
1. Use LQ_markov_mapping to map 𝑈𝑔,𝑡 , 𝐴22,𝑡 , 𝐶2,𝑡 , 𝑝𝑡,𝑡+1 , 𝑝𝑡,𝑡+2 into the 𝐴, 𝐵, 𝐶, 𝑅, 𝑄, 𝑊 matrices for each
of the 𝑛 Markov states.
2. Use the LQMarkov class to solve the resulting n-state Markov jump LQ problem.
To implement a simple example of the two-period model, we assume that 𝐺𝑡 follows an AR(1) process:
1
To do this, we set 𝑧𝑡 = [ ], and consequently:
𝐺𝑡
1 0 0
𝐴22 = [ ̄ ] , 𝐶2 = [ ] , 𝑈𝑔 = [0 1]
𝐺 𝜌 𝜎
We first solve the model with no penalty parameter on different issuance across maturities, i.e. 𝑐1 = 0.
We specify that the transition matrix for the Markov state is
0.9 0.1
Π=[ ]
0.1 0.9
Thus, each Markov state is persistent, and there is an equal chance of moving from one to the other.
# Model parameters
β, Gbar, ρ, σ, c1 = 0.95, 5, 0.8, 1, 0
p1, p2, p3, p4 = β, β**2 - 0.02, β, β**2 + 0.02
A1, B1, C1, R1, Q1, W1 = LQ_markov_mapping(A22, C_2, Ug, p1, p2, c1)
A2, B2, C2, R2, Q2, W2 = LQ_markov_mapping(A22, C_2, Ug, p3, p4, c1)
Π = np.array([[0.9, 0.1],
[0.1, 0.9]])
The above simulations show that when no penalty is imposed on different issuances across maturities, the government has
an incentive to take large “long-short” positions in debt of different maturities.
To prevent such outcomes, we set 𝑐1 = 0.01.
This penalty is big enough to motivate the government to issue positive quantities of both one- and two-period debt:
A1, B1, C1, R1, Q1, W1 = LQ_markov_mapping(A22, C_2, Ug, p1, p2, c1)
A2, B2, C2, R2, Q2, W2 = LQ_markov_mapping(A22, C_2, Ug, p3, p4, c1)
To map this into the Markov Jump LQ framework, we define state and control variables.
Let:
𝑏𝑡𝑡−1 𝑡
𝑏𝑡+1
⎡ 𝑏𝑡−1 ⎤ ⎡ 𝑏𝑡 ⎤
𝑏̄𝑡 = ⎢ 𝑡+1 ⎥ , 𝑢𝑡 = ⎢ 𝑡+2 ⎥
⎢ ⋮ ⎥ ⎢ ⋮ ⎥
𝑡−1 𝑡
𝑏
⎣ 𝑡+𝐻−1 ⎦ ⎣𝑏𝑡+𝐻 ⎦
Thus, 𝑏̄𝑡 is the endogenous state (debt issued last period) and 𝑢𝑡 is the control (debt issued today).
As before, we will also have the exogenous state 𝑧𝑡 , which determines government spending.
𝑏̄
𝑥𝑡 = [ 𝑡 ]
𝑧𝑡
We also define a vector 𝑝𝑡 that contains the time 𝑡 price of goods in period 𝑡 + 𝑗:
𝑝𝑡,𝑡+1
⎡𝑝 ⎤
𝑝𝑡 = ⎢ 𝑡,𝑡+2 ⎥
⎢ ⋮ ⎥
⎣𝑝𝑡,𝑡+𝐻 ⎦
𝑝𝑡,𝑡+1 1 0 0 ⋯ 0
⎡ 𝑝 ⎤ ⎡0 1 0 ⋯ 0⎤
⎢ 𝑡,𝑡+2 ⎥ = 𝑆𝑠 𝑝𝑡 where 𝑆𝑠 = ⎢ ⎥
⎢ ⋮ ⎥ ⎢⋮ ⋱ ⎥
⎣𝑝𝑡,𝑡+𝐻−1 ⎦ ⎣0 0 ⋯ 1 0⎦
𝑡−1
𝑏𝑡+1 0 1 0 ⋯ 0
⎡ 𝑏𝑡−1 ⎤ ⎡0 0 1 ⋯ 0⎤
⎢ 𝑡+2 ̄
⎥ = 𝑆𝑥 𝑏𝑡 where 𝑆𝑥 = ⎢ ⎥
⎢ ⋮ ⎥ ⎢⋮ ⋱ ⎥
𝑡−1
⎣𝑏𝑡+𝑇 −1 ⎦ ⎣0 0 ⋯ 0 1⎦
or
𝑇𝑡 = 𝑆𝑡 𝑥𝑡 − 𝑝𝑡′ 𝑢𝑡
Therefore
where
where to economize on notation we adopt the convention that for the linear state matrices 𝑅𝑡 ≡ 𝑅𝑠𝑡 , 𝑄𝑡 ≡ 𝑊𝑠𝑡 and so
on.
We’ll use this convention for the linear state matrices 𝐴, 𝐵, 𝑊 and so on below.
Because the payoff function also includes the penalty parameter for rescheduling, we have:
𝐻−1
𝑇𝑡2 + ∑ 𝑐2 (𝑏𝑡+𝑗
𝑡−1 𝑡
− 𝑏𝑡+𝑗+1 )2 = 𝑇𝑡2 + 𝑐2 (𝑏̄𝑡 − 𝑢𝑡 )′ (𝑏̄𝑡 − 𝑢𝑡 )
𝑗=0
Because the complete state is 𝑥𝑡 and not 𝑏̄𝑡 , we rewrite this as:
where 𝑆𝑐 = [𝐼 0]
Multiplying this out gives:
Therefore, with the cost term, we must amend our 𝑅, 𝑄, 𝑊 matrices as follows:
𝑅𝑡𝑐 = 𝑅𝑡 + 𝑐2 𝑆𝑐′ 𝑆𝑐
𝑄𝑐𝑡 = 𝑄𝑡 + 𝑐2 𝐼
𝑊𝑡𝑐 = 𝑊𝑡 − 𝑐2 𝑆𝑐
To finish mapping into the Markov jump LQ setup, we need to construct the law of motion for the full state.
This is simpler than in the previous setup, as we now have 𝑏̄𝑡+1 = 𝑢𝑡 .
Therefore:
𝑏̄𝑡+1
𝑥𝑡+1 ≡ [ ] = 𝐴𝑡 𝑥𝑡 + 𝐵𝑢𝑡 + 𝐶𝑡 𝑤𝑡+1
𝑧𝑡+1
where
0 0 𝐼 0
𝐴𝑡 = [ ], 𝐵 = [ ], 𝐶=[ ]
0 𝐴22,𝑡 0 𝐶2,𝑡
We can define a function that maps the primitives of the model with restructuring into the matrices required by the
LQMarkov class:
"""
Function which takes A22, C2, T, p_t, c and returns the
required matrices for the LQMarkov model: A, B, C, R, Q, W
Note, p_t should be a T by 1 matrix
(continues on next page)
# Create Sx, tSx, Ss, S_t matrices (tSx stands for \tilde S_x)
Ss = np.hstack((np.eye(T-1), np.zeros((T-1, 1))))
Sx = np.hstack((np.zeros((T-1, 1)), np.eye(T-1)))
tSx = np.zeros((1, T))
tSx[0, 0] = 1
# Create A, B, C matrices
A_T = np.hstack((np.zeros((T, T)), np.zeros((T, nz))))
A_B = np.hstack((np.zeros((nz, T)), A22))
A = np.vstack((A_T, A_B))
As an example let 𝐻 = 3.
Assume that there are two Markov states, one with a flatter yield curve, the other with a steeper yield curve.
In state 1, prices are:
1 1 1
𝑝𝑡,𝑡+1 = 0.9695 , 𝑝𝑡,𝑡+2 = 0.902 , 𝑝𝑡,𝑡+3 = 0.8369
We specify the same transition matrix and 𝐺𝑡 process that we used earlier.
A1, B1, C1, R1, Q1, W1 = LQ_markov_mapping_restruct(A22, C_2, Ug, H, p1, c2)
A2, B2, C2, R2, Q2, W2 = LQ_markov_mapping_restruct(A22, C_2, Ug, H, p2, c2)
fig, ax = plt.subplots()
ax.plot((u[0, :] / (u[0, :] + u[1, :] + u[2, :])))
ax.set_title('One-period debt issuance share')
ax.set_xlabel('Time')
plt.show()
ELEVEN
11.1 Overview
This lecture presents another application of Markov jump linear quadratic dynamic programming and constitutes a sequel
to an earlier lecture.
We again use a method introduced in lecture Markov Jump LQ dynamic programming to implement some ideas of [Barro,
1999] and [Barro and McCleary, 2003]) that extend the classic [Barro, 1979] model of tax smoothing.
[Barro, 1979] is about a government that borrows and lends in order to help it minimize an intertemporal measure of
distortions caused by taxes.
Technically, [Barro, 1979] looks a lot like a consumption-smoothing model.
Our generalization will also look like a souped-up consumption-smoothing model.
In this lecture, we describe a tax-smoothing problem of a government that faces roll-over risk.
In addition to what’s in Anaconda, this lecture deploys the quantecon library:
import quantecon as qe
import numpy as np
import matplotlib.pyplot as plt
Let 𝑇𝑡 denote tax collections, 𝛽 a discount factor, 𝑏𝑡,𝑡+1 time 𝑡 + 1 goods that the government promises to pay at 𝑡, 𝐺𝑡
𝑡
government purchases, 𝑝𝑡+1 the number of time 𝑡 goods received per time 𝑡 + 1 goods promised.
The stochastic process of government expenditures is exogenous.
The government’s problem is to choose a plan for borrowing and tax collections {𝑏𝑡+1 , 𝑇𝑡 }∞
𝑡=0 to minimize
∞
𝐸0 ∑ 𝛽 𝑡 𝑇𝑡2
𝑡=0
221
Advanced Quantitative Economics with Python
𝐺𝑡 = 𝑈𝑔,𝑡 𝑧𝑡
in Markov state 2.
Consequently, in the second Markov state, the government is unable to borrow, and the budget constraint becomes 𝑇𝑡 =
𝐺𝑡 + 𝑏𝑡−1,𝑡 .
However, if this is the only adjustment we make in our linear-quadratic model, the government will not set 𝑏𝑡,𝑡+1 = 0,
which is the outcome we want to express roll-over risk in period 𝑡.
Instead, the government would have an incentive to set 𝑏𝑡,𝑡+1 to a large negative number in state 2 – it would accumulate
large amounts of assets to bring into period 𝑡 + 1 because that is cheap
• Riccati equations will tell us this
Thus, we must represent “roll-over risk” some other way.
To force the government to set 𝑏𝑡,𝑡+1 = 0, we can instead extend the model to have four Markov states:
1. Good today, good yesterday
2. Good today, bad yesterday
3. Bad today, good yesterday
4. Bad today, bad yesterday
where good is a state in which effectively the government can issue debt and bad is a state in which effectively the
government can’t issue debt.
We’ll explain what effectively means shortly.
We now set
𝑡
𝑝𝑡+1 =𝛽
in all states.
In addition – and this is important because it defines what we mean by effectively – we put a large penalty on the 𝑏𝑡−1,𝑡
element of the state vector in states 2 and 4.
This will prevent the government from wishing to issue any debt in states 3 or 4 because it would experience a large
penalty from doing so in the next period.
The transition matrix for this formulation is:
0.95 0 0.05 0
⎡0.95 0 0.05 0 ⎤
Π=⎢ ⎥
⎢ 0 0.9 0 0.1⎥
⎣ 0 0.9 0 0.1⎦
This transition matrix ensures that the Markov state cannot move, for example, from state 3 to state 1.
Because state 3 is “bad today”, the next period cannot have “good yesterday”.
# Model parameters
β, Gbar, ρ, σ = 0.95, 5, 0.8, 1
# LQ framework matrices
A_t = np.zeros((1, 3))
A_b = np.hstack((np.zeros((2, 1)), A22))
A = np.vstack((A_t, A_b))
B = np.zeros((3, 1))
B[0, 0] = 1
R = S.T @ S
M = np.array([[-β]])
(continues on next page)
Using the same process for 𝐺𝑡 as in this lecture, we shall simulate our model with roll-over risk.
𝑡
When 𝑝𝑡+1 = 𝛽 government debt fluctuates around zero.
The spikes in the tax collection series indicate periods when the government is unable to access financial markets:
• positive spikes occur when debt is positive and the government must urgently raise tax revenues now
Negative spikes occur when the government has positive asset holdings.
An inability to use financial markets in the next period means that the government uses those assets to lower taxation
today.
x0 = np.array([[0, 1, 25]])
T = 300
x, u, w, state = lqm.compute_sequence(x0, ts_length=T)
# Calculate taxation each period from the budget constraint and the Markov state
tax = np.zeros([T, 1])
for i in range(T):
tax[i, :] = S @ x[:, i] + M @ u[:, i]
We can adjust parameters so that, rather than debt fluctuating around zero, the government is a debtor in every period
that it can borrow.
𝑡
To accomplish this, we simply raise 𝑝𝑡+1 to 𝛽 + 0.02 = 0.97.
M = np.array([[-β - 0.02]])
Q = M.T @ M
W = M.T @ S
# Calculate taxation each period from the budget constraint and the
# Markov state
tax = np.zeros([T, 1])
for i in range(T):
tax[i, :] = S @ x[:, i] + M @ u[:, i]
With a lower interest rate, the government has an incentive to increase debt over time.
However, with “roll-over risk”, debt is recurrently reset to zero and tax collections spike up.
In this model, high costs of a “sudden stop” make the government wary about letting its debt get too high.
TWELVE
In addition to what’s in Anaconda, this lecture will need the following libraries:
12.1 Overview
import sys
import numpy as np
import matplotlib.pyplot as plt
(continues on next page)
227
Advanced Quantitative Economics with Python
We begin by outlining the key assumptions regarding technology, households and the government sector.
12.2.1 Technology
12.2.2 Households
Consider a representative household who chooses a path {ℓ𝑡 , 𝑐𝑡 } for labor and consumption to maximize
1 ∞
−𝔼 ∑ 𝛽 𝑡 [(𝑐𝑡 − 𝑏𝑡 )2 + ℓ𝑡2 ] (12.1)
2 𝑡=0
Here
• 𝛽 is a discount factor in (0, 1).
• 𝑝𝑡0 is a scaled Arrow-Debreu price at time 0 of history contingent goods at time 𝑡 + 𝑗.
• 𝑏𝑡 is a stochastic preference parameter.
• 𝑑𝑡 is an endowment process.
• 𝜏𝑡 is a flat tax rate on labor income.
𝛽 𝑡 𝑝𝑡0
𝜋𝑡0 (𝑥𝑡 )
Thus, our scaled Arrow-Debreu price is the ordinary Arrow-Debreu price multiplied by the discount factor 𝛽 𝑡 and divided
by an appropriate probability.
The budget constraint (12.2) requires that the present value of consumption be restricted to equal the present value of
endowments, labor income and coupon payments on bond holdings.
12.2.3 Government
The government imposes a linear tax on labor income, fully committing to a stochastic path of tax rates at time zero.
The government also issues state-contingent debt.
Given government tax and borrowing plans, we can construct a competitive equilibrium with distorting government taxes.
Among all such competitive equilibria, the Ramsey plan is the one that maximizes the welfare of the representative
consumer.
Endowments, government expenditure, the preference shock process 𝑏𝑡 , and promised coupon payments on initial gov-
ernment debt 𝑠𝑡 are all exogenous, and given by
• 𝑑𝑡 = 𝑆𝑑 𝑥𝑡
• 𝑔𝑡 = 𝑆𝑔 𝑥𝑡
• 𝑏𝑡 = 𝑆𝑏 𝑥𝑡
• 𝑠𝑡 = 𝑆𝑠 𝑥𝑡
The matrices 𝑆𝑑 , 𝑆𝑔 , 𝑆𝑏 , 𝑆𝑠 are primitives and {𝑥𝑡 } is an exogenous stochastic process taking values in ℝ𝑘 .
We consider two specifications for {𝑥𝑡 }.
1. Discrete case: {𝑥𝑡 } is a discrete state Markov chain with transition matrix 𝑃 .
2. VAR case: {𝑥𝑡 } obeys 𝑥𝑡+1 = 𝐴𝑥𝑡 + 𝐶𝑤𝑡+1 where {𝑤𝑡 } is independent zero-mean Gaussian with identify
covariance matrix.
12.2.5 Feasibility
𝑐𝑡 + 𝑔𝑡 = 𝑑𝑡 + ℓ𝑡 (12.3)
Where 𝑝𝑡0 is again a scaled Arrow-Debreu price, the time zero government budget constraint is
∞
𝔼 ∑ 𝛽 𝑡 𝑝𝑡0 (𝑠𝑡 + 𝑔𝑡 − 𝜏𝑡 ℓ𝑡 ) = 0 (12.4)
𝑡=0
12.2.7 Equilibrium
An equilibrium is a feasible allocation {ℓ𝑡 , 𝑐𝑡 }, a sequence of prices {𝑝𝑡0 }, and a tax system {𝜏𝑡 } such that
1. The allocation {ℓ𝑡 , 𝑐𝑡 } is optimal for the household given {𝑝𝑡0 } and {𝜏𝑡 }.
2. The government’s budget constraint (12.4) is satisfied.
The Ramsey problem is to choose the equilibrium {ℓ𝑡 , 𝑐𝑡 , 𝜏𝑡 , 𝑝𝑡0 } that maximizes the household’s welfare.
If {ℓ𝑡 , 𝑐𝑡 , 𝜏𝑡 , 𝑝𝑡0 } solves the Ramsey problem, then {𝜏𝑡 } is called the Ramsey plan.
The solution procedure we adopt is
1. Use the first-order conditions from the household problem to pin down prices and allocations given {𝜏𝑡 }.
2. Use these expressions to rewrite the government budget constraint (12.4) in terms of exogenous variables and
allocations.
3. Maximize the household’s objective function (12.1) subject to the constraint constructed in step 2 and the feasibility
constraint (12.3).
The solution to this maximization problem pins down all quantities of interest.
12.2.8 Solution
Step one is to obtain the first-conditions for the household’s problem, taking taxes and prices as given.
Letting 𝜇 be the Lagrange multiplier on (12.2), the first-order conditions are 𝑝𝑡0 = (𝑐𝑡 − 𝑏𝑡 )/𝜇 and ℓ𝑡 = (𝑐𝑡 − 𝑏𝑡 )(1 − 𝜏𝑡 ).
Rearranging and normalizing at 𝜇 = 𝑏0 − 𝑐0 , we can write these conditions as
𝑏𝑡 − 𝑐 𝑡 ℓ𝑡
𝑝𝑡0 = and 𝜏𝑡 = 1 − (12.5)
𝑏0 − 𝑐 0 𝑏𝑡 − 𝑐 𝑡
The Ramsey problem now amounts to maximizing (12.1) subject to (12.6) and (12.3).
The associated Lagrangian is
∞
1
ℒ = 𝔼 ∑ 𝛽 𝑡 {− [(𝑐𝑡 − 𝑏𝑡 )2 + ℓ𝑡2 ] + 𝜆 [(𝑏𝑡 − 𝑐𝑡 )(ℓ𝑡 − 𝑠𝑡 − 𝑔𝑡 ) − ℓ𝑡2 ] + 𝜇𝑡 [𝑑𝑡 + ℓ𝑡 − 𝑐𝑡 − 𝑔𝑡 ]} (12.7)
𝑡=0
2
and
ℓ𝑡 − 𝜆[(𝑏𝑡 − 𝑐𝑡 ) − 2ℓ𝑡 ] = 𝜇𝑡
Combining these last two equalities with (12.3) and working through the algebra, one can show that
where
• 𝜈 ∶= 𝜆/(1 + 2𝜆)
• ℓ𝑡̄ ∶= (𝑏𝑡 − 𝑑𝑡 + 𝑔𝑡 )/2
• 𝑐𝑡̄ ∶= (𝑏𝑡 + 𝑑𝑡 − 𝑔𝑡 )/2
• 𝑚𝑡 ∶= (𝑏𝑡 − 𝑑𝑡 − 𝑠𝑡 )/2
Apart from 𝜈, all of these quantities are expressed in terms of exogenous variables.
To solve for 𝜈, we can use the government’s budget constraint again.
The term inside the brackets in (12.6) is (𝑏𝑡 − 𝑐𝑡 )(𝑠𝑡 + 𝑔𝑡 ) − (𝑏𝑡 − 𝑐𝑡 )ℓ𝑡 + ℓ𝑡2 .
Using (12.8), the definitions above and the fact that ℓ ̄ = 𝑏 − 𝑐,̄ this term can be rewritten as
𝑏0 + 𝑎0 (𝜈 2 − 𝜈) = 0
for 𝜈.
Provided that 4𝑏0 < 𝑎0 , there is a unique solution 𝜈 ∈ (0, 1/2), and a unique corresponding 𝜆 > 0.
Let’s work out how to compute mathematical expectations in (12.10).
For the first one, the random variable (𝑏𝑡 − 𝑐𝑡̄ )(𝑔𝑡 + 𝑠𝑡 ) inside the summation can be expressed as
1 ′
𝑥 (𝑆 − 𝑆𝑑 + 𝑆𝑔 )′ (𝑆𝑔 + 𝑆𝑠 )𝑥𝑡
2 𝑡 𝑏
For the second expectation in (12.10), the random variable 2𝑚2𝑡 can be written as
1 ′
𝑥 (𝑆 − 𝑆𝑑 − 𝑆𝑠 )′ (𝑆𝑏 − 𝑆𝑑 − 𝑆𝑠 )𝑥𝑡
2 𝑡 𝑏
It follows that both objects of interest are special cases of the expression
∞
𝑞(𝑥0 ) = 𝔼 ∑ 𝛽 𝑡 𝑥′𝑡 𝐻𝑥𝑡 (12.11)
𝑡=0
Next, suppose that {𝑥𝑡 } is the discrete Markov process described above.
Suppose further that each 𝑥𝑡 takes values in the state space {𝑥1 , … , 𝑥𝑁 } ⊂ ℝ𝑘 .
Let ℎ ∶ ℝ𝑘 → ℝ be a given function, and suppose that we wish to evaluate
∞
𝑞(𝑥0 ) = 𝔼 ∑ 𝛽 𝑡 ℎ(𝑥𝑡 ) given 𝑥0 = 𝑥𝑗
𝑡=0
Here
• 𝑃 𝑡 is the 𝑡-th power of the transition matrix 𝑃 .
• ℎ is, with some abuse of notation, the vector (ℎ(𝑥1 ), … , ℎ(𝑥𝑁 )).
• (𝑃 𝑡 ℎ)[𝑗] indicates the 𝑗-th element of 𝑃 𝑡 ℎ.
It can be shown that (12.12) is in fact equal to the 𝑗-th element of the vector (𝐼 − 𝛽𝑃 )−1 ℎ.
This last fact is applied in the calculations below.
We are interested in tracking several other variables besides the ones described above.
To prepare the way for this, we define
𝑡
𝑏𝑡+𝑗 − 𝑐𝑡+𝑗
𝑝𝑡+𝑗 =
𝑏𝑡 − 𝑐 𝑡
as the scaled Arrow-Debreu time 𝑡 price of a history contingent claim on one unit of consumption at time 𝑡 + 𝑗.
These are prices that would prevail at time 𝑡 if markets were reopened at time 𝑡.
These prices are constituents of the present value of government obligations outstanding at time 𝑡, which can be expressed
as
∞
𝐵𝑡 ∶= 𝔼𝑡 ∑ 𝛽 𝑗 𝑝𝑡+𝑗
𝑡
(𝜏𝑡+𝑗 ℓ𝑡+𝑗 − 𝑔𝑡+𝑗 ) (12.13)
𝑗=0
Using our expression for prices and the Ramsey plan, we can also write 𝐵𝑡 as
∞ 2
(𝑏𝑡+𝑗 − 𝑐𝑡+𝑗 )(ℓ𝑡+𝑗 − 𝑔𝑡+𝑗 ) − ℓ𝑡+𝑗
𝐵𝑡 = 𝔼𝑡 ∑ 𝛽 𝑗
𝑗=0
𝑏𝑡 − 𝑐 𝑡
and
𝑡
𝐵𝑡 = (𝜏𝑡 ℓ𝑡 − 𝑔𝑡 ) + 𝛽𝐸𝑡 𝑝𝑡+1 𝐵𝑡+1 (12.14)
Define
𝑅𝑡−1 ∶= 𝔼𝑡 𝛽 𝑗 𝑝𝑡+1
𝑡
(12.15)
12.2.12 A Martingale
where 𝐸𝑡̃ is the conditional mathematical expectation taken with respect to a one-step transition density that has been
formed by multiplying the original transition density with the likelihood ratio
𝑡
𝑝𝑡+1
𝑚𝑡𝑡+1 = 𝑡
𝐸𝑡 𝑝𝑡+1
which asserts that {𝜋𝑡+1 } is a martingale difference sequence under the distorted probability measure, and that {Π𝑡 } is
a martingale under the distorted probability measure.
In the tax-smoothing model of Robert Barro [Barro, 1979], government debt is a random walk.
In the current model, government debt {𝐵𝑡 } is not a random walk, but the excess payoff {Π𝑡 } on it is.
12.3 Implementation
Parameters
===========
T: int
Length of the simulation
Returns
========
path: a namedtuple of type 'Path', containing
g - Govt spending
d - Endowment
b - Utility shift parameter
s - Coupon payment on existing debt
c - Consumption
l - Labor
p - Price
τ - Tax rate
rvn - Revenue
B - Govt debt
R - Risk-free gross return
π - One-period risk-free interest rate
Π - Cumulative rate of return, adjusted
ξ - Adjustment factor for Π
"""
# Simplify names
β, Sg, Sd, Sb, Ss = econ.β, econ.Sg, econ.Sd, econ.Sb, econ.Ss
if econ.discrete:
P, x_vals = econ.proc
else:
A, C = econ.proc
return path
def gen_fig_1(path):
"""
The parameter is the path namedtuple returned by compute_paths(). See
the docstring of that function for details.
"""
T = len(path.c)
# Prepare axes
num_rows, num_cols = 2, 2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(14, 10))
plt.subplots_adjust(hspace=0.4)
for i in range(num_rows):
for j in range(num_cols):
axes[i, j].grid()
axes[i, j].set_xlabel('Time')
bbox = (0., 1.02, 1., .102)
legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}
p_args = {'lw': 2, 'alpha': 0.7}
plt.show()
def gen_fig_2(path):
"""
The parameter is the path namedtuple returned by compute_paths(). See
the docstring of that function for details.
"""
T = len(path.c)
# Prepare axes
num_rows, num_cols = 2, 1
fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 10))
plt.subplots_adjust(hspace=0.5)
plt.show()
The function var_quadratic_sum imported from quadsums is for computing the value of (12.11) when the ex-
ogenous process {𝑥𝑡 } is of the VAR type described above.
Below the definition of the function, you will see definitions of two namedtuple objects, Economy and Path.
The first is used to collect all the parameters and primitives of a given LQ economy, while the second collects output of
the computations.
In Python, a namedtuple is a popular data type from the collections module of the standard library that replicates
the functionality of a tuple, but also allows you to assign a name to each tuple element.
These elements can then be references via dotted attribute notation — see for example the use of path in the functions
gen_fig_1() and gen_fig_2().
The benefits of using namedtuples:
• Keeps content organized by meaning.
• Helps reduce the number of global variables.
Other than that, our code is long but relatively straightforward.
12.4 Examples
# == Parameters == #
β = 1 / 1.05
ρ, mg = .7, .35
A = eye(2)
A[0, :] = ρ, mg * (1-ρ)
C = np.zeros((2, 1))
C[0, 0] = np.sqrt(1 - ρ**2) * mg / 10
Sg = np.array((1, 0)).reshape(1, 2)
Sd = np.array((0, 0)).reshape(1, 2)
Sb = np.array((0, 2.135)).reshape(1, 2)
Ss = np.array((0, 0)).reshape(1, 2)
T = 50
path = compute_paths(T, economy)
gen_fig_1(path)
gen_fig_2(path)
Our second example adopts a discrete Markov specification for the exogenous process
# == Parameters == #
β = 1 / 1.05
P = np.array([[0.8, 0.2, 0.0],
[0.0, 0.5, 0.5],
[0.0, 0.0, 1.0]])
Sg = np.array((1, 0, 0, 0, 0)).reshape(1, 5)
Sd = np.array((0, 1, 0, 0, 0)).reshape(1, 5)
Sb = np.array((0, 0, 1, 0, 0)).reshape(1, 5)
Ss = np.array((0, 0, 0, 1, 0)).reshape(1, 5)
T = 15
path = compute_paths(T, economy)
gen_fig_1(path)
↪extract a single element from your array before performing this operation.␣
gen_fig_2(path)
12.5 Exercises
Exercise 12.5.1
Modify the VAR example given above, setting
# == Parameters == #
β = 1 / 1.05
ρ, mg = .95, .35
A = np.array([[0, 0, 0, ρ, mg*(1-ρ)],
[1, 0, 0, 0, 0],
[0, 1, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 1]])
C = np.zeros((5, 1))
C[0, 0] = np.sqrt(1 - ρ**2) * mg / 8
Sg = np.array((1, 0, 0, 0, 0)).reshape(1, 5)
Sd = np.array((0, 0, 0, 0, 0)).reshape(1, 5)
# Chosen st. (Sc + Sg) * x0 = 1
Sb = np.array((0, 0, 0, 0, 2.135)).reshape(1, 5)
Ss = np.array((0, 0, 0, 0, 0)).reshape(1, 5)
T = 50
path = compute_paths(T, economy)
gen_fig_1(path)
gen_fig_2(path)
249
CHAPTER
THIRTEEN
In addition to what’s in Anaconda, this lecture will need the following libraries:
13.1 Overview
This lecture computes versions of Arellano’s [Arellano, 2008] model of sovereign default.
The model describes interactions among default risk, output, and an equilibrium interest rate that includes a premium for
endogenous default risk.
The decision maker is a government of a small open economy that borrows from risk-neutral foreign creditors.
The foreign lenders must be compensated for default risk.
The government borrows and lends abroad in order to smooth the consumption of its citizens.
The government repays its debt only if it wants to, but declining to pay has adverse consequences.
The interest rate on government debt adjusts in response to the state-dependent default probability chosen by government.
The model yields outcomes that help interpret sovereign default experiences, including
• countercyclical interest rates on sovereign debt
• countercyclical trade balances
• high volatility of consumption relative to output
Notably, long recessions caused by bad draws in the income process increase the government’s incentive to default.
This can lead to
• spikes in interest rates
• temporary losses of access to international credit markets
• large drops in output, consumption, and welfare
• large capital outflows during recessions
Such dynamics are consistent with experiences of many countries.
Let’s start with some imports:
251
Advanced Quantitative Economics with Python
13.2 Structure
A small open economy is endowed with an exogenous stochastically fluctuating potential output stream {𝑦𝑡 }.
Potential output is realized only in periods in which the government honors its sovereign debt.
The output good can be traded or consumed.
The sequence {𝑦𝑡 } is described by a Markov process with stochastic density kernel 𝑝(𝑦, 𝑦′ ).
Households within the country are identical and rank stochastic consumption streams according to
∞
𝔼 ∑ 𝛽 𝑡 𝑢(𝑐𝑡 ) (13.1)
𝑡=0
Here
• 0 < 𝛽 < 1 is a time discount factor
• 𝑢 is an increasing and strictly concave utility function
Consumption sequences enjoyed by households are affected by the government’s decision to borrow or lend internationally.
The government is benevolent in the sense that its aim is to maximize (13.1).
The government is the only domestic actor with access to foreign credit.
Because household are averse to consumption fluctuations, the government will try to smooth consumption by borrowing
from (and lending to) foreign creditors.
The only credit instrument available to the government is a one-period bond traded in international credit markets.
The bond market has the following features
• The bond matures in one period and is not state contingent.
• A purchase of a bond with face value 𝐵′ is a claim to 𝐵′ units of the consumption good next period.
• To purchase 𝐵′ next period costs 𝑞𝐵′ now, or, what is equivalent.
• For selling −𝐵′ units of next period goods the seller earns −𝑞𝐵′ of today’s goods.
– If 𝐵′ < 0, then −𝑞𝐵′ units of the good are received in the current period, for a promise to repay −𝐵′ units
next period.
– There is an equilibrium price function 𝑞(𝐵′ , 𝑦) that makes 𝑞 depend on both 𝐵′ and 𝑦.
Earnings on the government portfolio are distributed (or, if negative, taxed) lump sum to households.
When the government is not excluded from financial markets, the one-period national budget constraint is
Here and below, a prime denotes a next period value or a claim maturing next period.
To rule out Ponzi schemes, we also require that 𝐵 ≥ −𝑍 in every period.
• 𝑍 is chosen to be sufficiently large that the constraint never binds in equilibrium.
Foreign creditors
• are risk neutral
• know the domestic output stochastic process {𝑦𝑡 } and observe 𝑦𝑡 , 𝑦𝑡−1 , … , at time 𝑡
• can borrow or lend without limit in an international credit market at a constant international interest rate 𝑟
• receive full payment if the government chooses to pay
• receive zero if the government defaults on its one-period debt due
When a government is expected to default next period with probability 𝛿, the expected value of a promise to pay one unit
of consumption next period is 1 − 𝛿.
Therefore, the discounted expected value of a promise to pay 𝐵 next period is
1−𝛿
𝑞= (13.3)
1+𝑟
Next we turn to how the government in effect chooses the default probability 𝛿.
While in a state of default, the economy regains access to foreign credit in each subsequent period with probability 𝜃.
13.3 Equilibrium
Informally, an equilibrium is a sequence of interest rates on its sovereign debt, a stochastic sequence of government default
decisions and an implied flow of household consumption such that
1. Consumption and assets satisfy the national budget constraint.
2. The government maximizes household utility taking into account
• the resource constraint
• the effect of its choices on the price of bonds
• consequences of defaulting now for future net output and future borrowing and lending opportunities
3. The interest rate on the government’s debt includes a risk-premium sufficient to make foreign creditors expect on
average to earn the constant risk-free international interest rate.
To express these ideas more precisely, consider first the choices of the government, which
1. enters a period with initial assets 𝐵, or what is the same thing, initial debt to be repaid now of −𝐵
2. observes current output 𝑦, and
3. chooses either
1. to default, or
2. to pay −𝐵 and set next period’s debt due to −𝐵′
In a recursive formulation,
• state variables for the government comprise the pair (𝐵, 𝑦)
• 𝑣(𝐵, 𝑦) is the optimum value of the government’s problem when at the beginning of a period it faces the choice of
whether to honor or default
• 𝑣𝑐 (𝐵, 𝑦) is the value of choosing to pay obligations falling due
• 𝑣𝑑 (𝑦) is the value of choosing to default
𝑣𝑑 (𝑦) does not depend on 𝐵 because, when access to credit is eventually regained, net foreign assets equal 0.
Expressed recursively, the value of defaulting is
𝑣𝑐 (𝐵, 𝑦) = max
′
{𝑢(𝑦 − 𝑞(𝐵′ , 𝑦)𝐵′ + 𝐵) + 𝛽 ∫ 𝑣(𝐵′ , 𝑦′ )𝑝(𝑦, 𝑦′ )𝑑𝑦′ }
𝐵 ≥−𝑍
Given zero profits for foreign creditors in equilibrium, we can combine (13.3) and (13.4) to pin down the bond price
function:
1 − 𝛿(𝐵′ , 𝑦)
𝑞(𝐵′ , 𝑦) = (13.5)
1+𝑟
An equilibrium is
• a pricing function 𝑞(𝐵′ , 𝑦),
• a triple of value functions (𝑣𝑐 (𝐵, 𝑦), 𝑣𝑑 (𝑦), 𝑣(𝐵, 𝑦)),
• a decision rule telling the government when to default and when to pay as a function of the state (𝐵, 𝑦), and
• an asset accumulation rule that, conditional on choosing not to default, maps (𝐵, 𝑦) into 𝐵′
such that
• The three Bellman equations for (𝑣𝑐 (𝐵, 𝑦), 𝑣𝑑 (𝑦), 𝑣(𝐵, 𝑦)) are satisfied
• Given the price function 𝑞(𝐵′ , 𝑦), the default decision rule and the asset accumulation decision rule attain the
optimal value function 𝑣(𝐵, 𝑦), and
• The price function 𝑞(𝐵′ , 𝑦) satisfies equation (13.5)
13.4 Computation
class Arellano_Economy:
" Stores data and creates primitives for the Arellano economy. "
def __init__(self,
B_grid_size= 251, # Grid size for bonds
B_grid_min=-0.45, # Smallest B value
B_grid_max=0.45, # Largest B value
y_grid_size=51, # Grid size for income
β=0.953, # Time discount parameter
γ=2.0, # Utility parameter
r=0.017, # Lending rate
ρ=0.945, # Persistence in the income process
η=0.025, # Standard deviation of the income process
θ=0.282, # Prob of re-entering financial markets
def_y_param=0.969): # Parameter governing income in default
# Save parameters
self.β, self.γ, self.r, = β, γ, r
self.ρ, self.η, self.θ = ρ, η, θ
self.y_grid_size = y_grid_size
self.B_grid_size = B_grid_size
self.B_grid = np.linspace(B_grid_min, B_grid_max, B_grid_size)
mc = qe.markov.tauchen(y_grid_size, ρ, η, 0, 3)
self.y_grid, self.P = np.exp(mc.state_values), mc.P
def params(self):
return self.β, self.γ, self.r, self.ρ, self.η, self.θ
def arrays(self):
return self.P, self.y_grid, self.B_grid, self.def_y, self.B0_idx
Notice how the class returns the data it stores as simple numerical values and arrays via the methods params and
arrays.
We will use this data in the Numba-jitted functions defined below.
Jitted functions prefer simple arguments, since type inference is easier.
Here is the utility function.
@njit
def u(c, γ):
return c**(1-γ)/(1-γ)
Here is a function to compute the bond price at each state, given 𝑣𝑐 and 𝑣𝑑 .
@njit
def compute_q(v_c, v_d, q, params, arrays):
"""
Compute the bond price function q(b, y) at each (b, y) pair.
# Unpack
β, γ, r, ρ, η, θ = params
P, y_grid, B_grid, def_y, B0_idx = arrays
@njit
def T_d(y_idx, v_c, v_d, params, arrays):
"""
The RHS of the Bellman equation when income is at index y_idx and
the country has chosen to default. Returns an update of v_d.
"""
# Unpack
β, γ, r, ρ, η, θ = params
P, y_grid, B_grid, def_y, B0_idx = arrays
current_utility = u(def_y[y_idx], γ)
v = np.maximum(v_c[B0_idx, :], v_d)
cont_value = np.sum((θ * v + (1 - θ) * v_d) * P[y_idx, :])
@njit
def T_c(B_idx, y_idx, v_c, v_d, q, params, arrays):
"""
The RHS of the Bellman equation when the country is not in a
defaulted state on their debt. Returns a value that corresponds to
v_c[B_idx, y_idx], as well as the optimal level of bond sales B'.
"""
# Unpack
β, γ, r, ρ, η, θ = params
P, y_grid, B_grid, def_y, B0_idx = arrays
B = B_grid[B_idx]
y = y_grid[y_idx]
Here is a fast function that calls these operators in the right sequence.
@njit(parallel=True)
def update_values_and_prices(v_c, v_d, # Current guess of value functions
B_star, q, # Arrays to be written to
params, arrays):
# Unpack
β, γ, r, ρ, η, θ = params
P, y_grid, B_grid, def_y, B0_idx = arrays
y_grid_size = len(y_grid)
B_grid_size = len(B_grid)
# Allocate memory
new_v_c = np.empty_like(v_c)
new_v_d = np.empty_like(v_d)
We can now write a function that will use the Arellano_Economy class and the functions defined above to compute
the solution to our model.
We do not need to JIT compile this function since it only consists of outer loops (and JIT compiling makes almost zero
difference).
In fact, one of the jobs of this function is to take an instance of Arellano_Economy, which is hard for the JIT
compiler to handle, and strip it down to more basic objects, which are then passed out to jitted functions.
# Allocate memory
q = np.empty_like(v_c)
B_star = np.empty_like(v_c, dtype=int)
current_iter = 0
dist = np.inf
while (current_iter < max_iter) and (dist > tol):
if current_iter % 100 == 0:
print(f"Entering iteration {current_iter}.")
Finally, we write a function that will allow us to simulate the economy once we have the policy functions
"""
# Unpack elements of the model
B0_idx = model.B0_idx
y_grid = model.y_grid
B_grid, y_grid, P = model.B_grid, model.y_grid, model.P
# Perform simulation
t = 0
while t < T:
# if in default:
if v_c[B_idx, y_idx] < v_d[y_idx] or in_default:
y_a_sim[t] = model.def_y[y_idx]
d_sim[t] = 1
Bp_idx = B0_idx
# Re-enter financial markets next period with prob θ
in_default = False if np.random.rand() < model.θ else True
else:
y_a_sim[t] = y_sim[t]
d_sim[t] = 0
Bp_idx = B_star[B_idx, y_idx]
13.5 Results
The grid used to compute this figure was relatively fine (y_grid_size, B_grid_size = 51, 251), which
explains the minor differences between this and Arrelano’s figure.
The figure shows that
• Higher levels of debt (larger −𝐵′ ) induce larger discounts on the face value, which correspond to higher interest
rates.
• Lower income also causes more discounting, as foreign creditors anticipate greater likelihood of default.
The next figure plots value functions and replicates the right hand panel of Figure 4 of [Arellano, 2008].
We can use the results of the computation to study the default probability 𝛿(𝐵′ , 𝑦) defined in (13.4).
The next plot shows these default probabilities over (𝐵′ , 𝑦) as a heat map.
As anticipated, the probability that the government chooses to default in the following period increases with indebtedness
and falls with income.
Next let’s run a time series simulation of {𝑦𝑡 }, {𝐵𝑡 } and 𝑞(𝐵𝑡+1 , 𝑦𝑡 ).
The grey vertical bars correspond to periods when the economy is excluded from financial markets because of a past
default.
One notable feature of the simulated data is the nonlinear response of interest rates.
Periods of relative stability are followed by sharp spikes in the discount rate on government debt.
13.6 Exercises
Exercise 13.6.1
To the extent that you can, replicate the figures shown above
• Use the parameter values listed as defaults in Arellano_Economy.
• The time series will of course vary depending on the shock draws.
ae = Arellano_Economy()
Entering iteration 0.
# Create "Y High" and "Y Low" values as 5% devs from mean
high, low = np.mean(y_grid) * 1.05, np.mean(y_grid) * .95
iy_high, iy_low = (np.searchsorted(y_grid, x) for x in (high, low))
# Create figure
fig, ax = plt.subplots(figsize=(10, 6.5))
hm = ax.pcolormesh(xx, yy, zz.T)
cax = fig.add_axes([.92, .1, .02, .8])
fig.colorbar(hm, cax=cax)
ax.axis([xx.min(), 0.05, yy.min(), yy.max()])
ax.set(xlabel="$B'$", ylabel="$y$", title="Probability of Default")
plt.show()
T = 250
np.random.seed(42)
y_sim, y_a_sim, B_sim, q_sim, d_sim = simulate(ae, T, v_c, v_d, q, B_star)
plt.show()
FOURTEEN
14.1 Overview
In this lecture, we review the paper Globalization and Synchronization of Innovation Cycles by Kiminori Matsuyama,
Laura Gardini and Iryna Sushko.
This model helps us understand several interesting stylized facts about the world economy.
One of these is synchronized business cycles across different countries.
Most existing models that generate synchronized business cycles do so by assumption, since they tie output in each country
to a common shock.
They also fail to explain certain features of the data, such as the fact that the degree of synchronization tends to increase
with trade ties.
By contrast, in the model we consider in this lecture, synchronization is both endogenous and increasing with the extent
of trade integration.
In particular, as trade costs fall and international competition increases, innovation incentives become aligned and coun-
tries synchronize their innovation cycles.
Let’s start with some imports:
import numpy as np
import matplotlib.pyplot as plt
from numba import jit
from ipywidgets import interact
14.1.1 Background
The model builds on work by Judd [Judd, 1985], Deneckner and Judd [Deneckere and Judd, 1992] and Helpman and
Krugman [Helpman and Krugman, 1985] by developing a two-country model with trade and innovation.
On the technical side, the paper introduces the concept of coupled oscillators to economic modeling.
As we will see, coupled oscillators arise endogenously within the model.
Below we review the model and replicate some of the results on synchronization of innovation across countries.
271
Advanced Quantitative Economics with Python
As discussed above, two countries produce and trade with each other.
In each country, firms innovate, producing new varieties of goods and, in doing so, receiving temporary monopoly power.
Imitators follow and, after one period of monopoly, what had previously been new varieties now enter competitive pro-
duction.
Firms have incentives to innovate and produce new goods when the mass of varieties of goods currently in production is
relatively low.
In addition, there are strategic complementarities in the timing of innovation.
Firms have incentives to innovate in the same period, so as to avoid competing with substitutes that are competitively
produced.
This leads to temporal clustering in innovations in each country.
After a burst of innovation, the mass of goods currently in production increases.
However, goods also become obsolete, so that not all survive from period to period.
This mechanism generates a cycle, where the mass of varieties increases through simultaneous innovation and then falls
through obsolescence.
14.2.2 Synchronization
In the absence of trade, the timing of innovation cycles in each country is decoupled.
This will be the case when trade costs are prohibitively high.
If trade costs fall, then goods produced in each country penetrate each other’s markets.
As illustrated below, this leads to synchronization of business cycles across the two countries.
14.3 Model
𝑜
Here 𝑋𝑘,𝑡 is a homogeneous input which can be produced from labor using a linear, one-for-one technology.
It is freely tradeable, competitively supplied, and homogeneous across countries.
By choosing the price of this good as numeraire and assuming both countries find it optimal to always produce the
homogeneous good, we can set 𝑤1,𝑡 = 𝑤2,𝑡 = 1.
The good 𝑋𝑘,𝑡 is a composite, built from many differentiated goods via
1
1− 1 1− 𝜎
𝑋𝑘,𝑡 𝜎 = ∫ [𝑥𝑘,𝑡 (𝜈)] 𝑑𝜈
Ω𝑡
Here 𝑥𝑘,𝑡 (𝜈) is the total amount of a differentiated good 𝜈 ∈ Ω𝑡 that is produced.
The parameter 𝜎 > 1 is the direct partial elasticity of substitution between a pair of varieties and Ω𝑡 is the set of varieties
available in period 𝑡.
We can split the varieties into those which are supplied competitively and those supplied monopolistically; that is, Ω𝑡 =
Ω𝑐𝑡 + Ω𝑚
𝑡 .
14.3.1 Prices
The price of a variety also depends on the origin, 𝑗, and destination, 𝑘, of the goods because shipping varieties between
countries incurs an iceberg trade cost 𝜏𝑗,𝑘 .
Thus the effective price in country 𝑘 of a variety 𝜈 produced in country 𝑗 becomes 𝑝𝑘,𝑡 (𝜈) = 𝜏𝑗,𝑘 𝑝𝑗,𝑡 (𝜈).
Using these expressions, we can derive the total demand for each variety, which is
where
𝜌𝑗,𝑘 𝐿𝑘
𝐴𝑗,𝑡 ∶= ∑ and 𝜌𝑗,𝑘 = (𝜏𝑗,𝑘 )1−𝜎 ≤ 1
𝑘
(𝑃𝑘,𝑡 )1−𝜎
It is assumed that 𝜏1,1 = 𝜏2,2 = 1 and 𝜏1,2 = 𝜏2,1 = 𝜏 for some 𝜏 > 1, so that
Monopolists will have the same marked-up price, so, for all 𝜈 ∈ Ω𝑚 ,
𝑚 𝜓 𝑚 𝑚 −𝜎
𝑝𝑗,𝑡 (𝜈) = 𝑝𝑗,𝑡 ∶= 1 and 𝐷𝑗,𝑡 = 𝑦𝑗,𝑡 ∶= 𝛼𝐴𝑗,𝑡 (𝑝𝑗,𝑡 )
1− 𝜎
Define
𝑐
𝑝𝑗,𝑡 𝑐
𝑦𝑗,𝑡 1 1−𝜎
𝜃 ∶= 𝑚 𝑚 = (1 − )
𝑝𝑗,𝑡 𝑦𝑗,𝑡 𝜎
Using the preceding definitions and some algebra, the price indices can now be rewritten as
1−𝜎 𝑚
𝑃𝑘,𝑡 𝑐
𝑁𝑗,𝑡
( ) = 𝑀𝑘,𝑡 + 𝜌𝑀𝑗,𝑡 where 𝑀𝑗,𝑡 ∶= 𝑁𝑗,𝑡 +
𝜓 𝜃
𝑐 𝑚
The symbols 𝑁𝑗,𝑡 and 𝑁𝑗,𝑡 will denote the measures of Ω𝑐 and Ω𝑚 respectively.
To introduce a new variety, a firm must hire 𝑓 units of labor per variety in each country.
Monopolist profits must be less than or equal to zero in expectation, so
𝑚 𝑚 𝑚 𝑚 𝑚 𝑚
𝑁𝑗,𝑡 ≥ 0, 𝜋𝑗,𝑡 ∶= (𝑝𝑗,𝑡 − 𝜓)𝑦𝑗,𝑡 −𝑓 ≤0 and 𝜋𝑗,𝑡 𝑁𝑗,𝑡 =0
𝑚 𝑐 1 𝛼𝐿𝑗 𝛼𝐿𝑘
𝑁𝑗,𝑡 = 𝜃(𝑀𝑗,𝑡 − 𝑁𝑗,𝑡 ) ≥ 0, [ + ]≤𝑓
𝜎 𝜃(𝑀𝑗,𝑡 + 𝜌𝑀𝑘,𝑡 ) 𝜃(𝑀𝑗,𝑡 + 𝑀𝑘,𝑡 /𝜌)
With 𝛿 as the exogenous probability of a variety becoming obsolete, the dynamic equation for the measure of firms
becomes
𝑐 𝑐 𝑚 𝑐 𝑐
𝑁𝑗,𝑡+1 = 𝛿(𝑁𝑗,𝑡 + 𝑁𝑗,𝑡 ) = 𝛿(𝑁𝑗,𝑡 + 𝜃(𝑀𝑗,𝑡 − 𝑁𝑗,𝑡 ))
Here
𝐷𝐿𝐿 ∶= {(𝑛1 , 𝑛2 ) ∈ ℝ2+ |𝑛𝑗 ≤ 𝑠𝑗 (𝜌)}
𝐷𝐻𝐻 ∶= {(𝑛1 , 𝑛2 ) ∈ ℝ2+ |𝑛𝑗 ≥ ℎ𝑗 (𝑛𝑘 )}
𝐷𝐻𝐿 ∶= {(𝑛1 , 𝑛2 ) ∈ ℝ2+ |𝑛1 ≥ 𝑠1 (𝜌) and 𝑛2 ≤ ℎ2 (𝑛1 )}
𝐷𝐿𝐻 ∶= {(𝑛1 , 𝑛2 ) ∈ ℝ2+ |𝑛1 ≤ ℎ1 (𝑛2 ) and 𝑛2 ≥ 𝑠2 (𝜌)}
while
𝑠1 − 𝜌𝑠2
𝑠1 (𝜌) = 1 − 𝑠2 (𝜌) = min { , 1}
1−𝜌
14.4 Simulation
@jit(nopython=True)
def _hj(j, nk, s1, s2, θ, δ, ρ):
"""
If we expand the implicit function for h_j(n_k) then we find that
it is quadratic. We know that h_j(n_k) > 0 so we can get its
value by using the quadratic form
"""
# Find out who's h we are evaluating
if j == 1:
sj = s1
sk = s2
else:
sj = s2
(continues on next page)
return root
@jit(nopython=True)
def DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
"Determine whether (n1, n2) is in the set DLL"
return (n1 <= s1_ρ) and (n2 <= s2_ρ)
@jit(nopython=True)
def DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
"Determine whether (n1, n2) is in the set DHH"
return (n1 >= _hj(1, n2, s1, s2, θ, δ, ρ)) and \
(n2 >= _hj(2, n1, s1, s2, θ, δ, ρ))
@jit(nopython=True)
def DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
"Determine whether (n1, n2) is in the set DHL"
return (n1 >= s1_ρ) and (n2 <= _hj(2, n1, s1, s2, θ, δ, ρ))
@jit(nopython=True)
def DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
"Determine whether (n1, n2) is in the set DLH"
return (n1 <= _hj(1, n2, s1, s2, θ, δ, ρ)) and (n2 >= s2_ρ)
@jit(nopython=True)
def one_step(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
"""
Takes a current value for (n_{1, t}, n_{2, t}) and returns the
values (n_{1, t+1}, n_{2, t+1}) according to the law of motion.
"""
# Depending on where we are, evaluate the right branch
if DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
n1_tp1 = δ * (θ * s1_ρ + (1 - θ) * n1)
n2_tp1 = δ * (θ * s2_ρ + (1 - θ) * n2)
elif DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
n1_tp1 = δ * n1
n2_tp1 = δ * n2
elif DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
n1_tp1 = δ * n1
n2_tp1 = δ * (θ * _hj(2, n1, s1, s2, θ, δ, ρ) + (1 - θ) * n2)
elif DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ):
n1_tp1 = δ * (θ * _hj(1, n2, s1, s2, θ, δ, ρ) + (1 - θ) * n1)
n2_tp1 = δ * n2
@jit(nopython=True)
@jit(nopython=True)
def _pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, maxiter, npers):
"""
Takes initial values and iterates forward to see whether
the histories eventually end up in sync.
If countries are symmetric then as soon as the two countries have the
same measure of firms then they will be synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm
Parameters
----------
n1_0 : scalar(Float)
Initial normalized measure of firms in country one
n2_0 : scalar(Float)
Initial normalized measure of firms in country two
maxiter : scalar(Int)
Maximum number of periods to simulate
npers : scalar(Int)
Number of periods we would like the countries to have the
same measure for
Returns
-------
synchronized : scalar(Bool)
Did the two economies end up synchronized
pers_2_sync : scalar(Int)
The number of periods required until they synchronized
"""
# Initialize the status of synchronization
synchronized = False
pers_2_sync = maxiter
iters = 0
# Initialize generator
n_gen = n_generator(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
@jit(nopython=True)
def _create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ,
maxiter, npers, npts):
# Create unit range with npts
synchronized, pers_2_sync = False, 0
unit_range = np.linspace(0.0, 1.0, npts)
return time_2_sync
class MSGSync:
"""
The paper "Globalization and Synchronization of Innovation Cycles" presents
a two-country model with endogenous innovation cycles. Combines elements
from Deneckere Judd (1985) and Helpman Krugman (1985) to allow for a
model with trade that has firms who can introduce new varieties into
the economy.
Parameters
----------
s1 : scalar(Float)
def _unpack_params(self):
return self.s1, self.s2, self.θ, self.δ, self.ρ
def _calc_s1_ρ(self):
# Unpack params
s1, s2, θ, δ, ρ = self._unpack_params()
# s_1(ρ) = min(val, 1)
val = (s1 - ρ * s2) / (1 - ρ)
return min(val, 1)
Parameters
----------
n1_0 : scalar(Float)
Initial normalized measure of firms in country one
n2_0 : scalar(Float)
Initial normalized measure of firms in country two
T : scalar(Int)
Number of periods to simulate
Returns
-------
n1 : Array(Float64, ndim=1)
A history of normalized measures of firms in country one
n2 : Array(Float64, ndim=1)
A history of normalized measures of firms in country two
"""
# Unpack parameters
s1, s2, θ, δ, ρ = self._unpack_params()
s1_ρ, s2_ρ = self.s1_ρ, self.s2_ρ
# Allocate space
n1 = np.empty(T)
n2 = np.empty(T)
# Store in arrays
n1[t] = n1_tp1
n2[t] = n2_tp1
return n1, n2
If countries are symmetric then as soon as the two countries have the
same measure of firms then they will be synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm
Parameters
----------
n1_0 : scalar(Float)
Initial normalized measure of firms in country one
n2_0 : scalar(Float)
Initial normalized measure of firms in country two
maxiter : scalar(Int)
Maximum number of periods to simulate
npers : scalar(Int)
Number of periods we would like the countries to have the
same measure for
Returns
-------
synchronized : scalar(Bool)
Did the two economies end up synchronized
pers_2_sync : scalar(Int)
The number of periods required until they synchronized
"""
# Unpack parameters
s1, s2, θ, δ, ρ = self._unpack_params()
s1_ρ, s2_ρ = self.s1_ρ, self.s2_ρ
return ab
We write a short function below that exploits the preceding code and plots two time series.
Each time series gives the dynamics for the two countries.
The time series share parameters but differ in their initial condition.
Here’s the function
ax.legend()
ax.set(title=title, ylim=(0.15, 0.8))
return ax
# Create figure
fig, ax = plt.subplots(2, 1, figsize=(10, 8))
fig.tight_layout()
plt.show()
In the first case, innovation in the two countries does not synchronize.
In the second case, different initial conditions are chosen, and the cycles become synchronized.
Next, let’s study the initial conditions that lead to synchronized cycles more systematically.
We generate time series from a large collection of different initial conditions and mark those conditions with different
colors according to whether synchronization occurs or not.
The next display shows exactly this for four different parameterizations (one for each subfigure).
Dark colors indicate synchronization, while light colors indicate failure to synchronize.
As you can see, larger values of 𝜌 translate to more synchronization.
You are asked to replicate this figure in the exercises.
In the solution to the exercises, you’ll also find a figure with sliders, allowing you to experiment with different parameters.
Here’s one snapshot from the interactive figure
14.5 Exercises
Exercise 14.5.1
Replicate the figure shown above by coloring initial conditions according to whether or not synchronization occurs from
those conditions.
return ab, cf
Additionally, instead of just seeing 4 plots at once, we might want to manually be able to change 𝜌 and see how it affects
the plot in real-time. Below we use an interactive plot to do this.
Note, interactive plotting requires the ipywidgets module to be installed and enabled.
Note: This interactive plot is disabled on this static webpage. In order to use this, we recommend to run this notebook
locally.
fig = interact(interact_attraction_basis,
ρ=(0.0, 1.0, 0.05),
maxiter=(50, 5000, 50),
npts=(25, 750, 25))
FIFTEEN
15.1 Overview
In 1937, Ronald Coase wrote a brilliant essay on the nature of the firm [Coase, 1937].
Coase was writing at a time when the Soviet Union was rising to become a significant industrial power.
At the same time, many free-market economies were afflicted by a severe and painful depression.
This contrast led to an intensive debate on the relative merits of decentralized, price-based allocation versus top-down
planning.
In the midst of this debate, Coase made an important observation: even in free-market economies, a great deal of top-
down planning does in fact take place.
This is because firms form an integral part of free-market economies and, within firms, allocation is by planning.
In other words, free-market economies blend both planning (within firms) and decentralized production coordinated by
prices.
The question Coase asked is this: if prices and free markets are so efficient, then why do firms even exist?
Couldn’t the associated within-firm planning be done more efficiently by the market?
We’ll use the following imports:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fminbound
On top of asking a deep and fascinating question, Coase also supplied an illuminating answer: firms exist because of
transaction costs.
Here’s one example of a transaction cost:
Suppose agent A is considering setting up a small business and needs a web developer to construct and help run an online
store.
She can use the labor of agent B, a web developer, by writing up a freelance contract for these tasks and agreeing on a
suitable price.
But contracts like this can be time-consuming and difficult to verify
• How will agent A be able to specify exactly what she wants, to the finest detail, when she herself isn’t sure how the
business will evolve?
289
Advanced Quantitative Economics with Python
• And what if she isn’t familiar with web technology? How can she specify all the relevant details?
• And, if things go badly, will failure to comply with the contract be verifiable in court?
In this situation, perhaps it will be easier to employ agent B under a simple labor contract.
The cost of this contract is far smaller because such contracts are simpler and more standard.
The basic agreement in a labor contract is: B will do what A asks him to do for the term of the contract, in return for a
given salary.
Making this agreement is much easier than trying to map every task out in advance in a contract that will hold up in a
court of law.
So agent A decides to hire agent B and a firm of nontrivial size appears, due to transaction costs.
15.1.2 A Trade-Off
15.1.3 Summary
15.2.1 Subcontracting
The subcontracting scheme by which tasks are allocated across firms is illustrated in the figure below
In this example,
• Firm 1 receives a contract to sell one unit of the completed good to a final buyer.
• Firm 1 then forms a contract with firm 2 to purchase the partially completed good at stage 𝑡1 , with the intention of
implementing the remaining 1 − 𝑡1 tasks in-house (i.e., processing from stage 𝑡1 to stage 1).
• Firm 2 repeats this procedure, forming a contract with firm 3 to purchase the good at stage 𝑡2 .
• Firm 3 decides to complete the chain, selecting 𝑡3 = 0.
At this point, production unfolds in the opposite direction (i.e., from upstream to downstream).
• Firm 3 completes processing stages from 𝑡3 = 0 up to 𝑡2 and transfers the good to firm 2.
• Firm 2 then processes from 𝑡2 up to 𝑡1 and transfers the good to firm 1,
• Firm 1 processes from 𝑡1 to 1 and delivers the completed good to the final buyer.
The length of the interval of stages (range of tasks) carried out by firm 𝑖 is denoted by ℓ𝑖 .
Each firm chooses only its upstream boundary, treating its downstream boundary as given.
The benefit of this formulation is that it implies a recursive structure for the decision problem for each firm.
In choosing how many processing stages to subcontract, each successive firm faces essentially the same decision problem
as the firm above it in the chain, with the only difference being that the decision space is a subinterval of the decision
space for the firm above.
We will exploit this recursive structure in our study of equilibrium.
15.2.2 Costs
15.3 Equilibrium
We assume that all firms are ex-ante identical and act as price takers.
As price takers, they face a price function 𝑝, which is a map from [0, 1] to ℝ+ , with 𝑝(𝑡) interpreted as the price of the
good at processing stage 𝑡.
There is a countable infinity of firms indexed by 𝑖 and no barriers to entry.
The cost of supplying the initial input (the good processed up to stage zero) is set to zero for simplicity.
Free entry and the infinite fringe of competitors rule out positive profits for incumbents, since any incumbent could be
replaced by a member of the competitive fringe filling the same role in the production chain.
Profits are never negative in equilibrium because firms can freely exit.
An equilibrium in this setting is an allocation of firms and a price function such that
1. all active firms in the chain make zero profits, including suppliers of raw materials
2. no firm in the production chain has an incentive to deviate, and
3. no inactive firms can enter and extract positive profits
In particular, 𝑡𝑖−1 is the downstream boundary of firm 𝑖 and 𝑡𝑖 is its upstream boundary.
As transaction costs are incurred only by the buyer, its profits are
1. 𝑝(0) = 0,
2. 𝜋𝑖 = 0 for all 𝑖, and
3. 𝑝(𝑠) − 𝑐(𝑠 − 𝑡) − 𝛿𝑝(𝑡) ≤ 0 for any pair 𝑠, 𝑡 with 0 ≤ 𝑠 ≤ 𝑡 ≤ 1.
The rationale behind these conditions was given in our informal definition of equilibrium above.
We have defined an equilibrium but does one exist? Is it unique? And, if so, how can we compute it?
To address these questions, we introduce the operator 𝑇 mapping a nonnegative function 𝑝 on [0, 1] to 𝑇 𝑝 via
𝑇 𝑝(𝑠) = min {𝑐(𝑠 − 𝑡) + 𝛿𝑝(𝑡)} for all 𝑠 ∈ [0, 1]. (15.3)
𝑡≤𝑠
By definition, 𝑡∗ (𝑠) is the cost-minimizing upstream boundary for a firm that is contracted to deliver the good at stage 𝑠
and faces the price function 𝑝∗ .
Since 𝑝∗ lies in 𝒫 and since 𝑐 is strictly convex, it follows that the right-hand side of (15.4) is continuous and strictly
convex in 𝑡.
Hence the minimizer 𝑡∗ (𝑠) exists and is uniquely defined.
We can use 𝑡∗ to construct an equilibrium allocation as follows:
Recall that firm 1 sells the completed good at stage 𝑠 = 1, its optimal upstream boundary is 𝑡∗ (1).
Hence firm 2’s optimal upstream boundary is 𝑡∗ (𝑡∗ (1)).
Continuing in this way produces the sequence {𝑡∗𝑖 } defined by
The sequence ends when a firm chooses to complete all remaining tasks.
We label this firm (and hence the number of firms in the chain) as
The task allocation corresponding to (15.5) is given by ℓ𝑖∗ ∶= 𝑡∗𝑖−1 − 𝑡∗𝑖 for all 𝑖.
In [Kikuchi et al., 2018] it is shown that
1. The value 𝑛∗ in (15.6) is well-defined and finite,
2. the allocation {ℓ𝑖∗ } is feasible, and
3. the price function 𝑝∗ and this allocation together forms an equilibrium for the production chain.
While the proofs are too long to repeat here, much of the insight can be obtained by observing that, as a fixed point of
𝑇 , the equilibrium price function must satisfy
From this equation, it is clear that so profits are zero for all incumbent firms.
We can develop some additional insights on the behavior of firms by examining marginal conditions associated with the
equilibrium.
As a first step, let ℓ∗ (𝑠) ∶= 𝑠 − 𝑡∗ (𝑠).
This is the cost-minimizing range of in-house tasks for a firm with downstream boundary 𝑠.
In [Kikuchi et al., 2018] it is shown that 𝑡∗ and ℓ∗ are increasing and continuous, while 𝑝∗ is continuously differentiable
at all 𝑠 ∈ (0, 1) with
Equation (15.8) follows from 𝑝∗ (𝑠) = min𝑡≤𝑠 {𝑐(𝑠 − 𝑡) + 𝛿𝑝∗ (𝑡)} and the envelope theorem for derivatives.
A related equation is the first order condition for 𝑝∗ (𝑠) = min𝑡≤𝑠 {𝑐(𝑠 − 𝑡) + 𝛿𝑝∗ (𝑡)}, the minimization problem for a
firm with upstream boundary 𝑠, which is
This condition matches the marginal condition expressed verbally by Coase that we stated above:
“A firm will tend to expand until the costs of organizing an extra transaction within the firm become equal
to the costs of carrying out the same transaction by means of an exchange on the open market…”
Combining (15.8) and (15.9) and evaluating at 𝑠 = 𝑡𝑖 , we see that active firms that are adjacent satisfy
𝛿 𝑐′ (ℓ𝑖+1
∗
) = 𝑐′ (ℓ𝑖∗ ) (15.10)
In other words, the marginal in-house cost per task at a given firm is equal to that of its upstream partner multiplied by
gross transaction cost.
This expression can be thought of as a Coase–Euler equation, which determines inter-firm efficiency by indicating how
two costly forms of coordination (markets and management) are jointly minimized in equilibrium.
15.5 Implementation
For most specifications of primitives, there is no closed-form solution for the equilibrium as far as we are aware.
However, we know that we can compute the equilibrium corresponding to a given transaction cost parameter 𝛿 and a cost
function 𝑐 by applying the results stated above.
In particular, we can
1. fix initial condition 𝑝 ∈ 𝒫,
2. iterate with 𝑇 until 𝑇 𝑛 𝑝 has converged to 𝑝∗ , and
3. recover firm choices via the choice function (15.3)
At each iterate, we will use continuous piecewise linear interpolation of functions.
To begin, here’s a class to store primitives and a grid:
class ProductionChain:
def __init__(self,
n=1000,
delta=1.05,
c=lambda t: np.exp(10 * t) - 1):
* pc is an instance of ProductionChain
* The initial condition is p = c
"""
delta, c, n, grid = pc.delta, pc.c, pc.n, pc.grid
p = c(grid) # Initial condition is c(s), as an array
new_p = np.empty_like(p)
error = tol + 1
i = 0
if i < max_iter:
print(f"Iteration converged in {i} steps")
else:
print(f"Warning: iteration hit upper bound {max_iter}")
(continues on next page)
The next function computes optimal choice of upstream boundary and range of task implemented for a firm face price
function p_function and with downstream boundary 𝑠.
"""
delta, c = pc.delta, pc.c
f = lambda t: delta * p_function(t) + c(s - t)
t_star = max(fminbound(f, -1, s), 0)
ell_star = s - t_star
return t_star, ell_star
The allocation of firms can be computed by recursively stepping through firms’ choices of their respective upstream
boundary, treating the previous firm’s upstream boundary as their own downstream boundary.
In doing so, we start with firm 1, who has downstream boundary 𝑠 = 1.
pc = ProductionChain()
p_star = compute_prices(pc)
fig, ax = plt.subplots()
ax.plot(pc.grid, p_star(pc.grid))
ax.set_xlim(0.0, 1.0)
ax.set_ylim(0.0)
for s in transaction_stages:
ax.axvline(x=s, c="0.5")
plt.show()
Here’s the function ℓ∗ , which shows how large a firm with downstream boundary 𝑠 chooses to be
ell_star = np.empty(pc.n)
for i, s in enumerate(pc.grid):
t, e = optimal_choices(pc, p_star, s)
ell_star[i] = e
fig, ax = plt.subplots()
ax.plot(pc.grid, ell_star, label=r"$\ell^*$")
ax.legend(fontsize=14)
plt.show()
15.6 Exercises
Exercise 15.6.1
The number of firms is endogenously determined by the primitives.
What do you think will happen in terms of the number of firms as 𝛿 increases? Why?
Check your intuition by computing the number of firms at delta in (1.01, 1.05, 1.1).
pc = ProductionChain(delta=delta)
p_star = compute_prices(pc)
transaction_stages = compute_stages(pc, p_star)
num_firms = len(transaction_stages)
print(f"When delta={delta} there are {num_firms} firms")
Exercise 15.6.2
The value added of firm 𝑖 is 𝑣𝑖 ∶= 𝑝∗ (𝑡𝑖−1 ) − 𝑝∗ (𝑡𝑖 ).
One of the interesting predictions of the model is that value added is increasing with downstreamness, as are several other
measures of firm size.
Can you give any intution?
Try to verify this phenomenon (value added increasing with downstreamness) using the code above.
pc = ProductionChain()
p_star = compute_prices(pc)
stages = compute_stages(pc, p_star)
va = []
fig, ax = plt.subplots()
ax.plot(va, label="value added by firm")
ax.set_xticks((5, 25))
ax.set_xticklabels(("downstream firms", "upstream firms"))
plt.show()
SIXTEEN
COMPOSITE SORTING
16.1 Overview
Optimal transport theory is studies how one (marginal) probabilty measure can be related to another (marginal) probability
measure in an ideal way.
The output of such a theory is a coupling of the two probability measures, i.e., a joint probabilty measure having those
two marginal probability measures.
This lecture describes how Job Boerma, Aleh Tsyvinski, Ruodo Wang, and Zhenyuan Zhang [Boerma et al., 2024] used
optimal transport theory to formulate and solve an equilibrium of a model in which wages and allocations of workers
across jobs adjust to match measures of different types with measures of different types of occupations.
Production technologies allow firms to affect shape costs of mismatch with the consequence that costs of mismatch can
be concave.
That means that it is possible that equilibrium there is neither positive assortive nor negative assorting matching, an
outcome that [Boerma et al., 2024] call composite assortive matching.
For example, in an equilibrium with composite matching, identical workers can sort into different occupations, some
positively and some negatively.
[Boerma et al., 2024] show how this can generate distinct distributions of labor earnings within and across occupations.
This lecture describes the [Boerma et al., 2024] model and presents Python code for computing equilibria.
The lecture applies the code to the [Boerma et al., 2024] model of labor markets.
As with an earlier QuantEcon lecture on optimal transport, a key tool will be linear programming.
16.2 Setup
𝑋 and 𝑌 are finite sets that represent two distinct types of people to be matched.
For each 𝑥 ∈ 𝑋, let a positive integer 𝑛𝑥 be the number of agents of type 𝑥.
Similarly, let a positive integer 𝑚𝑦 be the agents of agents of type 𝑦 ∈ 𝑌 .
We refer to these two measures as marginals.
We assume that
∑ 𝑛𝑥 = ∑ 𝑚𝑦 =∶ 𝑁
𝑥∈𝑋 𝑦∈𝑌
303
Advanced Quantitative Economics with Python
s.t. ∑ 𝜇𝑥𝑦 = 𝑛𝑥
𝑥∈𝑋
∑ 𝜇𝑥𝑦 = 𝑚𝑦
𝑦∈𝑌
Given our discreteness assumptions about 𝑋 and 𝑌 , the problem admits an integer solution 𝜇 ∈ ℤ𝑋×𝑌
+ , i.e., 𝜇𝑥𝑦 is a
non-negative integer for each 𝑥 ∈ 𝑋, 𝑦 ∈ 𝑌 .
We will study integer solutions.
Two points about restricting ourselves to integer solutions are worth mentioning:
• it is without loss of generality for computational purposes, since every problem with float marginals can be trans-
formed into an equivalent problem with integer marginals;
• although the mathematical structure that we present actually works for arbitrary real marginals, some of our Python
implementations would fail to work with float arithmetic.
We focus on a specific instance of an optimal transport problem:
We assume that 𝑋 and 𝑌 are finite subsets of ℝ and that the cost function satisfies 𝑐𝑥𝑦 = ℎ(|𝑥 − 𝑦|) for all 𝑥, 𝑦 ∈ ℝ, for
an ℎ ∶ ℝ+ → ℝ+ that is strictly concave and strictly increasing and grounded (i.e., ℎ(0) = 0).
Such an ℎ satisfies the following
Lemma. If ℎ ∶ ℝ+ → ℝ+ is strictly concave and grounded, then ℎ is strictly subadditive, i.e. for all 𝑥, 𝑦 ∈ ℝ+ , 0 < 𝑥 < 𝑦,
we have
Proof. For 𝛼 ∈ (0, 1) and 𝑥 > 0 we have, by strict concavity and groundedness, ℎ(𝛼𝑥) > 𝛼ℎ(𝑥)+(1−𝛼)ℎ(0) = 𝛼ℎ(𝑥).
𝑥
Now fix 𝑥, 𝑦 ∈ ℝ+ , 0 < 𝑥 < 𝑦, and let 𝛼 = 𝑥+𝑦 ; the previous observation gives ℎ(𝑥) = ℎ(𝛼(𝑥 + 𝑦)) > 𝛼ℎ(𝑥 + 𝑦) and
ℎ(𝑦) = ℎ((1 − 𝛼)(𝑥 + 𝑦)) > (1 − 𝛼)ℎ(𝑥 + 𝑦); summing these inequality delivers the result. □
In the following implementation we assume that the cost function is 𝑐𝑥𝑦 = |𝑥 − 𝑦|1/𝜁 for 𝜁 > 1, i.e. ℎ(𝑧) = 𝑧 1/𝜁 for
𝑧 ∈ ℝ+ .
Hence, our problem is
s.t. ∑ 𝜇𝑥𝑦 = 𝑛𝑥
𝑥∈𝑋
∑ 𝜇𝑥𝑦 = 𝑚𝑦
𝑦∈𝑌
import numpy as np
from scipy.optimize import linprog
from itertools import chain
import pandas as pd
from collections import namedtuple
(continues on next page)
The following Python class takes as inputs sets of types 𝑋, 𝑌 ⊂ ℝ, marginals 𝑛, 𝑚 with positive integer entries such that
∑𝑥∈𝑋 𝑛𝑥 = ∑𝑦∈𝑌 𝑚𝑦 and cost parameter 𝜁 > 1.
The cost function is stored as an |𝑋| × |𝑌 | matrix with (𝑥, 𝑦)-entry equal to |𝑥 − 𝑦|1/𝜁 , i.e., the cost of matching an agent
of type 𝑥 ∈ 𝑋 with an agent of type 𝑦 ∈ 𝑌 .
class ConcaveCostOT():
def __init__(self, X_types=None, Y_types=None, n_x =None, m_y=None, ζ=2):
# Sets of types
self.X_types, self.Y_types = X_types, Y_types
# Marginals
if X_types is not None and Y_types is not None:
non_empty_types = True
self.n_x = np.ones(len(X_types), dtype=int) if n_x is None else n_x
self.m_y = np.ones(len(Y_types), dtype=int) if m_y is None else m_y
else:
non_empty_types = False
self.n_x, self.m_y = n_x, m_y
Let’s consider a random instance with given numbers of types |𝑋| and |𝑌 | and a given number of agents.
First, we generate random types 𝑋 and 𝑌 .
Then we generate random quantities for each type so that there are 𝑁 agents for each side.
number_of_x_types = 20
number_of_y_types = 20
N_agents_per_side = 60
np.random.seed(1)
# generate types
X_types_example = np.random.choice(random_support,
(continues on next page)
ConcaveCostOT.assign_random_marginals = assign_random_marginals
We use 𝐹 (resp. 𝐺) to denote the cumulative distribution function associated to the measure 𝑛 (resp. 𝑚)
Thus, 𝐹 (𝑧) = ∑𝑥≤𝑧∶𝑛 𝑛𝑥 and 𝐺(𝑧) = ∑𝑦≤𝑧∶𝑚 𝑚𝑦 for 𝑧 ∈ ℝ.
𝑥 >0 𝑦 >0
plt.figure(figsize=figsize)
ConcaveCostOT.plot_marginals = plot_marginals
example_pb.plot_marginals()
We can verify that reassigning the minimum of such quantities to the pairs (𝑧, 𝑧) and (𝑥, 𝑦) improves upon the current
matching since
where the first inequality follows from triangle inequality and the fact that ℎ is increasing and the strict inequality from
strict subadditivity.
We can then repeat the operation for any other analogous pair of matches involving 𝑧, while improving the value, until
we have mass min{𝑛𝑧 , 𝑚𝑧 } on match (𝑧, 𝑧).
Viewing the matching 𝜇 as a measure on 𝑋 × 𝑌 with marginals 𝑛 and 𝑚, this property says that in any optimal 𝜇 we
have 𝜇𝑧𝑧 = 𝑛𝑧 ∧ 𝑚𝑧 for (𝑧, 𝑧) in the diagonal {(𝑥, 𝑦) ∈ 𝑋 × 𝑌 ∶ 𝑥 = 𝑦} of ℝ × ℝ.
The following method finds perfect pairs and returns the on-diagonal matchings as well as the residual off-diagonal
marginals.
def match_perfect_pairs(self):
m_y_off_diag = self.m_y.copy()
m_y_off_diag[perfect_pairs_y] -= Δ_q
ConcaveCostOT.match_perfect_pairs = match_perfect_pairs
On-diagonal matches: 15
Residual types in X: 14
Residual types in Y: 16
We can therefore create a new instance with the residual marginals that will feature no perfect pairs.
Later we shall add the on-diagonal matching to the solution of this new instance.
We refer to this instance as “off-diagonal” since the product measure of the residual marginals 𝑛 ⊗ 𝑚 feature zeros mass
on the diagonal of ℝ × ℝ.
In the rest of this section, we will focus on this instance.
We create a subclass to study the residual off-diagonal problem.
The subclass inherits the attributes and the modules from the original class.
We let 𝑍 ∶= 𝑋 ⊔ 𝑌 , where ⊔ denotes the union of disjoint sets. We will
• index types 𝑋 as {0, … , |𝑋| − 1} and types 𝑌 as {|𝑋|, … , |𝑋| + |𝑌 | − 1};
• store the cost function as a |𝑍| × |𝑍| matrix with entry (𝑧, 𝑧 ′ ) equal to 𝑐𝑥𝑦 if 𝑧 = 𝑥 ∈ 𝑋 and 𝑧 ′ = 𝑦 ∈ 𝑌 or
𝑧 = 𝑦 ∈ 𝑌 and 𝑧 ′ = 𝑥 ∈ 𝑋 or equal to +∞ if 𝑧 and 𝑧 ′ belong to the same side
– (the latter is just customary, since these “infinitely penalized” entries are actually never accessed in the im-
plementation);
• let 𝑞 be a vector of size |𝑍| whose 𝑧-th entry equals 𝑛𝑥 if type 𝑥 is the 𝑧-th smallest type in 𝑍 and −𝑚𝑦 if type 𝑦
is the 𝑧-th smallest type in 𝑍; hence 𝑞 encodes capacities of both sides on the (ascending) sorted set of types.
Finally, we add a method to flexibly add a pair (𝑖, 𝑗) with 𝑖 ∈ {0, … , |𝑋| − 1}, 𝑗 ∈ {|𝑋|, … , |𝑋| + |𝑌 | − 1} or
𝑗 ∈ {0, … , |𝑋| − 1}, 𝑖 ∈ {|𝑋|, … , |𝑋| + |𝑌 | − 1} to a matching matrix of size |𝑋| × |𝑌 |.
class OffDiagonal(ConcaveCostOT):
def __init__(self, X_types, Y_types, n_x, m_y, ζ):
super().__init__(X_types, Y_types, n_x, m_y, ζ)
# Types (unsorted)
self.types_list = np.concatenate((X_types,Y_types))
# upper-right block
self.cost_z_z[:len(self.X_types), len(self.X_types):] = self.cost_x_y
# lower-left block
self.cost_z_z[len(self.X_types):, :len(self.X_types)] = self.cost_x_y.T
## Distributions of types
# sorted types and index identifier for each z in support
self.type_z = np.argsort(self.types_list)
self.support_z = self.types_list[self.type_z]
We add a function that returns an instance of the off-diagonal subclass as well as the on-diagonal matching and the indices
of the residual off-diagonal types.
These indices will come handy for adding the off-diagonal matching matrix to the diagonal matching matrix we just found,
since the former will have a smaller size if there are perfect pairs in the original problem.
def generate_offD_onD_matching(self):
# Match perfect pairs and compute on-diagonal matching
n_x_off_diag, m_y_off_diag , matching_diag = self.match_perfect_pairs()
ConcaveCostOT.generate_offD_onD_matching = generate_offD_onD_matching
example_off_diag, _ = example_pb.generate_offD_onD_matching()
Let’s plot the residual marginals to verify visually that there are no overlappings between types from distinct sides in the
off-diagonal instance.
|𝑥 − 𝑦′ | + |𝑥′ − 𝑦| = |𝑥 − 𝑦| + |𝑥′ − 𝑦′ |
|𝑥−𝑦|+|𝑥′ −𝑦|
Letting 𝛼 ∶= |𝑥−𝑦′ |−|𝑥′ −𝑦| ∈ (0, 1), we have |𝑥−𝑦| = 𝛼|𝑥−𝑦′ |+(1−𝛼)|𝑥′ −𝑦| and |𝑥′ −𝑦′ | = (1−𝛼)|𝑥−𝑦′ |+𝛼|𝑥′ −𝑦|.
Hence, by strict concavity of ℎ,
ℎ(|𝑥 − 𝑦|) + ℎ(|𝑥′ − 𝑦′ |) < 𝛼ℎ(|𝑥 − 𝑦′ |) + (1 − 𝛼)ℎ(|𝑥′ − 𝑦|) + (1 − 𝛼)ℎ(|𝑥 − 𝑦′ |) + 𝛼ℎ(|𝑥′ − 𝑦|) = ℎ(|𝑥 − 𝑦′ |) + ℎ(|𝑥′ − 𝑦|).
Therefore, as in the first case, we can strictly improve the cost among 𝑥, 𝑦, 𝑥′ , 𝑦′ by uncrossing the pairs.
Finally, it remains to argue that in both cases uncrossing operations do not increase the number of intersections with other
matched pairs.
It can indeed be shown on a case-by-case basis that, in both of the above cases, for any other matched pair (𝑥″ , 𝑦″ ) the
number of intersections between pairs (𝑥, 𝑦), (𝑥′ , 𝑦′ ) and the pair (𝑥″ , 𝑦″ ) (i.e., after uncrossing) is not larger than the
number of intersections between pairs (𝑥, 𝑦′ ), (𝑥′ , 𝑦) and the pair (𝑥″ , 𝑦″ ) (i.e., before uncrossing), hence the uncrossing
operations above reduce the number of intersections.
We conclude that if a matching features intersecting pairs, it can be modified via a sequence of uncrossing operations
into a matching without intersecting pairs while improving on the value.
(Layering) Recall that there are 2𝑁 individual agents, each agent 𝑖 having type 𝑧𝑖 ∈ 𝑋 ⊔ 𝑌 .
When we introduce the off diagonal matching, to stress that the types sets are disjoint now.
To simplify our explanation of this property, assume for now that each agent has its own distinct type (i.e., |𝑋| = |𝑌 | = 𝑁
and 𝑛 = 𝑚 = 1𝑁 ), in which case the optimal transport problem is also referred to as assignment problem.
Let’s index agents according to their types:
Suppose that agents 𝑖 of type 𝑧𝑖 and 𝑗 of type 𝑧𝑗 , with 𝑧𝑖 < 𝑧𝑗 , are matched in a particular optimal solution.
Then there is an equal number of agents from each side in {𝑖 + 1, … , 𝑗 − 1}, if this set is not empty.
Indeed, if this were not the case, then some agent 𝑘 ∈ {𝑖 + 1, 𝑗 − 1} would be matched with some agent ℓ with
ℓ ∉ {𝑖, … , 𝑗}, i.e., there would be types
with matches (𝑧𝑖 , 𝑧𝑗 ) and (𝑧𝑘 , 𝑧ℓ ), violating the no intersecting pairs property.
We conclude that we can define a binary relation on [𝑁 ] such that 𝑖 ∼ 𝑗 if there is an equal number of agents of each
side in {𝑖, 𝑖 + 1, … , 𝑗} (or if this set is empty).
This is an equivalence relation, so we can find associated equivalence classes that we call layers.
By the reasoning above, in an optimal solution all pairs 𝑖, 𝑗 (of opposite sides) which are matched belong to the same
layer, hence we can solve the assignment problem associated to each layer and then add up the solutions.
In terms of distributions, 𝑖 and 𝑗, of types 𝑥 ∈ 𝑋 and 𝑦 ∈ 𝑌 respectively, belong to the same layer (i.e., 𝑥 ∼ 𝑦) if and
only if 𝐹 (𝑦−) − 𝐹 (𝑥) = 𝐺(𝑦−) − 𝐺(𝑥).
If 𝐹 and 𝐺 were continuous, then 𝐹 (𝑦) − 𝐹 (𝑥) = 𝐺(𝑦) − 𝐺(𝑥) ⟺ 𝐹 (𝑥) − 𝐺(𝑥) = 𝐹 (𝑦) − 𝐺(𝑦).
This suggests that the following quantity plays an important role:
# Plot H(z)
plt.figure(figsize=figsize)
plt.axhline(0, color='black', linewidth=1)
OffDiagonal.plot_H_z = plot_H_z
example_off_diag.plot_H_z()
Moreover, each layer 𝐿ℓ contains an even number of types 𝑁ℓ ∈ 2ℕ, which are alternating, i.e., ordering them as
𝑧1 < 𝑧2 ⋯ < 𝑧𝑁ℓ −1 < 𝑧𝑁ℓ all odd (or even, respectively) indexed types belong to the same side.
The following method finds the layers associated with distributions 𝐹 and 𝐺.
Again, types in 𝑋 are indexed with {0, … , |𝑋| − 1} and types in 𝑌 with {|𝑋|, … , |𝑋| + |𝑌 | − 1}.
Using these indices (instead of the types themselves) to represent the layers allows keeping track of sides types in each
layer, without adding an additional bit of information that would identify the side of the first type in the layer, which,
because a layer is alternating, would then allow identifying sides of all types in the layer.
In addition, using indices will let us extract the cost function within a layer from the cost function 𝑐𝑧𝑧′ computed offline.
def find_layers(self):
# Compute H(z) on the joint support
H_z = np.concatenate([[0], np.cumsum(self.q_z)])
# Compute layers
# the following |H(R)|x|Z| matrix has entry (z,l) equal to 1 iff type z belongs␣
↪to layer l
OffDiagonal.find_layers = find_layers
[array([23, 10]), array([27, 3, 23, 10]), array([16, 2, 21, 3, 25, 8, 23, 12]),
↪ array([16, 2, 21, 3, 25, 12]), array([22, 0, 16, 2, 21, 3, 18, 12]),␣
↪array([15, 0, 16, 2, 14, 5, 21, 3, 18, 9]), array([20, 0, 16, 2, 14, 5,␣
↪21, 3, 19, 11, 24, 1, 18, 9]), array([ 2, 26, 5, 21, 3, 19, 4, 18]),␣
↪array([ 2, 26, 7, 21, 3, 19, 4, 17, 6, 18]), array([13, 26, 7, 21, 3, 19, ␣
↪6, 18]), array([ 6, 18]), array([ 6, 28]), array([ 6, 29])]
plt.figure(figsize=figsize)
# Plot H(z)
step = np.concatenate(([self.support_z.min() - .05 * self.support_z.ptp()],
self.support_z,
[self.support_z.max() + .05 * self.support_z.ptp()]))
height = np.concatenate((H_z, [0]))
plt.step(step, height, where='post', color='black', label='CDF', zorder=1)
# Plot layers
colors = cm.viridis(np.linspace(0, 1, len(layers)))
for ell, layer in enumerate(layers):
plt.vlines(self.types_list[layer], layers_height[ell] ,
layers_height[ell] + layers_mass[ell],
color=colors[ell], linewidth=2)
plt.scatter(self.types_list[layer],
np.ones(len(layer)) * layers_height[ell]
+.5 * layers_mass[ell],
color=colors[ell], s=50)
plt.axhline(layers_height[ell], color=colors[ell],
linestyle=':', linewidth=1.5, zorder=0)
OffDiagonal.plot_layers = plot_layers
example_off_diag.plot_layers()
which is alternating.
The problem within a layer is unitary.
Hence we can solve the problem with unit masses and later rescale the solution by the layer’s mass 𝑀ℓ .
Let us select a layer from the example above (we pick the one with maximum number of types) and plot the types on the
real line
plt.figure(figsize=figsize)
ConcaveCostOT.plot_layer_types = plot_layer_types
example_off_diag.plot_layer_types(layer_example,
layers_mass_example[layer_id_example])
Given the structure of a layer and the no intersecting pairs property, the optimal matching and value of the layer can be
found recursively.
Indeed, if in certain optimal matching 1 and 𝑗 ∈ [𝑁ℓ ], 𝑗 − 1 odd, are paired, then there is no matching between agents
in [2, 𝑗 − 1] and those in [𝑗 + 1, 𝑁ℓ ] (if both are non empty, i.e., 𝑗 is not 2 or 𝑁ℓ ).
Hence in such optimal solution agents in [2, 𝑗 − 1] are matched among themselves.
Since [𝑧2 , 𝑧𝑗−1 ] (as well as [𝑧𝑗+1 , 𝑧𝑁ℓ ]) is alternating, we can reason recursively.
Let 𝑉𝑖𝑗 be the optimal value of matching agents in [𝑖, 𝑗] with 𝑖, 𝑗 ∈ [𝑁ℓ ], 𝑗 − 𝑖 ∈ {1, 3, … , 𝑁ℓ − 1}.
Suppose that we computed the value 𝑉𝑖𝑗 for all 𝑖, 𝑗 ∈ [𝑁ℓ ] with 𝑖 − 𝑗 ∈ {1, 3, … , 𝑡 − 2} for some odd natural number 𝑡.
Then, for 𝑖, 𝑗 ∈ [𝑁ℓ ] with 𝑖 − 𝑗 = 𝑡 we have
The following method takes as input the layer types indices and computes the value function as a matrix
[𝑉𝑖𝑗 ]𝑖∈[𝑁ℓ +1],𝑗∈[𝑁ℓ ] .
In order to distinguish entries that are relevant for our computations from those that are never accessed, we initialize this
matrix as full of NaN values.
def solve_bellman_eqs(self,layer):
# Recover cost function within the layer
cost_i_j = self.cost_z_z[layer[:,None],layer[None,:]]
t = 1
while t < len(layer):
# Select agents i in [n_L-t] (with potential partners j's in [t,n_L])
i_t = np.arange(len(layer)-t)
return V_i_j
OffDiagonal.solve_bellman_eqs = solve_bellman_eqs
Having computed the value function, we can proceed to compute the optimal matching as the policy that attains the value
function that solves the Bellman equation (policy evaluation).
We start from agent 1 and match it with the 𝑘 that achieves the minimum in the equation associated with 𝑉1,2𝑁ℓ .
Then we store segments [2, 𝑘 − 1] and [𝑘 + 1, 2𝑁ℓ ] (if not empty).
In general, given a segment [𝑖, 𝑗], we match 𝑖 with 𝑘 that achieves the minimum in the equation associated with 𝑉𝑖𝑗 and
store the segments [𝑖, 𝑘 − 1] and [𝑘 + 1, 𝑗] (if not empty).
The algorithm proceeds until there are no segments left.
while segments_to_process:
# Pick i, first agent of the segment
# and potential partners i+1,i+3,..., in the segment
segment = segments_to_process[0]
i_0 = segment[0]
potential_matches = np.arange(i_0, segment[-1], 2) + 1
return matching
OffDiagonal.find_layer_matching = find_layer_matching
Lets apply this method our example to find the matching within the layer and then rescale it by 𝑀ℓ .
Note that the unscaled value equals 𝑉1,𝑁ℓ .
matching_layer = example_off_diag.find_layer_matching(V_i_j,layer_example)
print(f"Value of the layer (unscaled): {(matching_layer * example_off_diag.cost_x_y).
↪sum()}")
plt.show()
ConcaveCostOT.plot_layer_matching = plot_layer_matching
example_off_diag.plot_layer_matching(layer_example, matching_layer)
We now present two key results in the context of OT with concave type costs.
We refer [Boerma et al., 2024] and [Delon et al., 2011] for proofs.
Consider the problem faced within a layer, i.e., types from 𝑌 ⊔ 𝑋
for 𝑖, 𝑗 ∈ [𝑁ℓ ], 𝑗 − 𝑖 odd, with boundary conditions 𝑉𝑖+1,𝑖 = 0 for 𝑖 ∈ [0, 𝑁ℓ ] and 𝑉𝑖+2,𝑖−1 = −𝑐𝑖,𝑖+1 for 𝑖 ∈ [𝑁ℓ − 1] .
The following method uses these equations to compute the value function that is stored as a matrix [𝑉𝑖𝑗 ]𝑖∈[𝑁ℓ +1],𝑗∈[𝑁ℓ +1] .
def solve_bellman_eqs_DSS(self,layer):
# Recover cost function within the layer
cost_i_j = self.cost_z_z[layer[:,None],layer[None,:]]
t = 1
while t < len(layer):
# Select agents i in [n_l-t] and potential partner j=i+t for each i
i_t = np.arange(len(layer)-t)
j_t = i_t + t +1
return V_i_j
OffDiagonal.solve_bellman_eqs_DSS = solve_bellman_eqs_DSS
Let’s apply the algorithm to our example and compare outcomes with those attained with the Bellman equations above.
V_i_j_DSS = example_off_diag.solve_bellman_eqs_DSS(layer_example)
print('##########################')
print(f"Difference with previous Bellman equations: \
{(V_i_j_DSS[:,1:] - V_i_j)[V_i_j >= 0].sum()}")
We can actually compute the optimal matching within the layer simultaneously with computing the value function, rather
than sequentially.
The key idea is that, if at some step of the computation of the values the left branch of the minimum above achieves the
minimum, say 𝑉𝑖𝑗 = 𝑐𝑖𝑗 + 𝑉𝑖+1,𝑗−1 , then (𝑖, 𝑗) are optimally matched on [𝑖, 𝑗] and by the theorem above we get that a
matching on [𝑖 + 1, 𝑗 − 1] which achieves 𝑉𝑖+1,𝑗−1 belongs to an optimal matching on the whole layer (since it is covered
by the arc (𝑖, 𝑗) in [𝑖, 𝑗]).
We can therefore proceed as follows
We initialize an empty matching and a list with all the agents in the layer (representing the agents which are not matched
yet).
Then whenever the left branch of the minimum is achieved for some (𝑖, 𝑗) in the computation of 𝑉 , we take the collections
of agents 𝑘1 , … , 𝑘𝑀 in [𝑖 + 1, 𝑗 − 1] (in ascending order, i.e. with 𝑧𝑘𝑝 < 𝑧𝑘𝑝+1 ) that are not matched yet (if any) and
add to the matching the pairs (𝑘1 , 𝑘2 ), (𝑘3 , 𝑘4 ), … , (𝑘𝑀−1 , 𝑘𝑀 ).
Thus, we match each unmatched agent 𝑘𝑝 in [𝑖 + 1, 𝑗 − 1] with the closest unmatched right neighbour 𝑘𝑝+1 (starting from
𝑘1 ).
Intuitively, if 𝑘𝑝 were optimally matched with some 𝑘𝑞 in [𝑖 + 1, 𝑗 − 1] and not with 𝑘𝑝+1 , then 𝑘𝑝+1 would have already
been hidden by the match (𝑘𝑝 , 𝑘𝑞 ) from some previous computation (because |𝑘𝑝 − 𝑘𝑞 | < |𝑖 − 𝑗|) and it would therefore
be matched.
Finally, if the process above leaves some umatched agents, we proceed by matching each of these agent with the closest
unmatched right neighbour, starting again from the leftmost of these collection.
To gain understanding, note that this situation happens when the left branch is achieved only for pairs 𝑖, 𝑗 with |𝑖 − 𝑗| = 1,
which leads to the optimal matching (1, 2), (2, 3), … , (𝑛ℓ − 1, 𝑛ℓ ).
def find_layer_matching_DSS(self,layer):
# Recover cost function within the layer
cost_i_j = self.cost_z_z[layer[:,None],layer[None,:]]
t = 1
while t < len(layer):
# Compute optimal value for pairs with |i-j| = t
i_t = np.arange(len(layer)-t)
j_t = i_t + t + 1
# Select each i for which left branch achieves minimum in the V_{i,i+t}
left_branch_achieved = i_t[left_branch == V_i_j[i_t, j_t]]
# Update matching
for i in left_branch_achieved:
# for each agent k in [i+1,i+t-1]
for k in np.arange(i+1,i+t)[unmatched[range(i+1,i+t)]]:
# if k is unmatched
if unmatched[k] == True:
# find unmatched right neighbour
j_k = np.arange(k+1,len(layer))[unmatched[k+1:]][0]
# add pair to matching
self.add_pair_to_matching(layer[[k, j_k]], matching)
# remove pair from unmatched agents list
unmatched[[k, j_k]] = False
return matching
OffDiagonal.find_layer_matching_DSS = find_layer_matching_DSS
matching_layer_DSS = example_off_diag.find_layer_matching_DSS(layer_example)
print(f" Value of layer with DSS recursive equations \
{(matching_layer_DSS * example_off_diag.cost_x_y).sum()}")
print(f" Value of layer with Bellman equations \
{(matching_layer * example_off_diag.cost_x_y).sum()}")
example_off_diag.plot_layer_matching(layer_example, matching_layer_DSS)
The following method assembles our components in order to solve the primal problem.
First, if matches are perfect pairs, we store the on-diagonal matching and create an off-diagonal instance with the residual
marginals.
Then we compute the set of layers of the residual distributions.
Finally, we solve each layer and put together matchings within each layer with the on-diagonal matchings.
We then return the full matching, the off-diagonal matching, and the off-diagonal instance.
def solve_primal_pb(self):
# Compute on-diagonal matching, create new instance with resitual types
off_diagoff_diagonal, match_tuple = self.generate_offD_onD_matching()
nonzero_id_x, nonzero_id_y, matching_diag = match_tuple
# Compute layers
(continues on next page)
ConcaveCostOT.solve_primal_pb = solve_primal_pb
def solve_primal_DSS(self):
# Compute on-diagonal matching, create new instance with resitual types
off_diagoff_diagonal, match_tuple = self.generate_offD_onD_matching()
nonzero_id_x, nonzero_id_y, matching_diag = match_tuple
# Find layers
layers, layers_mass, _, _ = off_diagoff_diagonal.find_layers()
ConcaveCostOT.solve_primal_DSS = solve_primal_DSS
DSS_tuple = example_pb.solve_primal_DSS()
matching_DSS, matching_off_diag_DSS, off_diagoff_diagonal_DSS = DSS_tuple
By drawing semicircles joining the matched agents (with distinct types), we can visualize the off-diagonal matching.
In the following figure, widths and colors of semicirles indicate relative numbers of agents that are “transported” along
an arc.
# Add labels
for i, x in enumerate(self.X_types):
ax.annotate(f'$x_{{{i }}}$', (x, 0), textcoords="offset points",
xytext=(0, -15), ha='center', color='blue', fontsize=12)
for j, y in enumerate(self.Y_types):
ax.annotate(f'$y_{{{j }}}$', (y, 0), textcoords="offset points",
xytext=(0, -15), ha='center', color='red', fontsize=12)
count = matching_off_diag[matched_types]
colors = plt.cm.Greys(np.linspace(0.5, 1.5, count.max() + 1))
max_height = 0
for iter in range(len(count)):
width = abs(matched_types_x[iter] - matched_types_y[iter])
center = (matched_types_x[iter] + matched_types_y[iter]) / 2
height = width
max_height = max(max_height, height)
semicircle = patches.Arc((center, 0), width, height,
theta1=0, theta2=180,
color=colors[count[iter]],
lw=count[iter] * (2.2 / count.max()))
ax.add_patch(semicircle)
step = np.concatenate(([self.support_z.min()
- .02 * self.support_z.ptp()],
self.support_z,
[self.support_z.max()
+ .02 * self.support_z.ptp()]))
# Set the y-limit to keep H_z and maximum circle size in the plot
ax.set_ylim(np.min(H_z) - H_z.ptp() *.01,
np.maximum(np.max(H_z), max_height / 2) + H_z.ptp() *.01)
plt.show()
ConcaveCostOT.plot_matching = plot_matching
off_diagoff_diagonal.plot_matching(matching_off_diag,
tit