From 8bf3b2c77cb5c92d859dcc3eae3566c49743740c Mon Sep 17 00:00:00 2001 From: thomassargent30 Date: Sat, 26 Sep 2026 15:35:26 -0400 Subject: [PATCH] Add three lectures and revise auction and fault-tree lectures New lectures: * lucas_prescott_investment: industry equilibrium under uncertainty following Lucas and Prescott (1971), with the planning-problem route to equilibrium, the invariant distribution, and exercises. * optimal_growth_uncertainty: Brock and Mirman (1972) with irreversible investment and Tobin's q following Sargent (1980), including a corrected argument for differentiability of the planner's value function. * lq_mean_field_games: linear quadratic mean field games following Alvarez and Argente (2026), with quantitative implementations of their two examples -- an industry equilibrium with capital accumulation, and multiproduct price setting with Kimball demand -- and eight exercises. Revisions: * house_auction, two_auctions: correct the VCG computations, replace the greedy assignment with an exact one, and add exercises; move the list of papers into the bibliography. * hoist_failure: explain the roles of grid range and resolution, compare the computed quantiles with the published fault-tree study, and add exercises. * rational_expectations: point to the two new lectures. Co-Authored-By: Claude Opus 5 (1M context) --- lectures/_static/quant-econ.bib | 235 ++++ lectures/_toc.yml | 3 + lectures/hoist_failure.md | 535 ++++++-- lectures/house_auction.md | 593 ++++++--- lectures/lq_mean_field_games.md | 1644 ++++++++++++++++++++++++ lectures/lucas_prescott_investment.md | 989 ++++++++++++++ lectures/optimal_growth_uncertainty.md | 842 ++++++++++++ lectures/rational_expectations.md | 34 + lectures/two_auctions.md | 331 ++++- 9 files changed, 4907 insertions(+), 299 deletions(-) create mode 100644 lectures/lq_mean_field_games.md create mode 100644 lectures/lucas_prescott_investment.md create mode 100644 lectures/optimal_growth_uncertainty.md diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index aa6536d31..0a51f5991 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -5445,3 +5445,238 @@ @article{Geweke1982 pages = {304--313}, year = {1982} } + +@article{Blackwell1965, + author = {Blackwell, David}, + title = {Discounted Dynamic Programming}, + journal = {Annals of Mathematical Statistics}, + volume = {36}, + number = {1}, + pages = {226--235}, + year = {1965} +} + +@article{HansenSargent1980, + author = {Hansen, Lars Peter and Sargent, Thomas J.}, + title = {Formulating and Estimating Dynamic Linear Rational Expectations + Models}, + journal = {Journal of Economic Dynamics and Control}, + volume = {2}, + pages = {7--46}, + year = {1980} +} + +@article{BrockMirman1972, + author = {Brock, William A. and Mirman, Leonard J.}, + title = {Optimal Economic Growth and Uncertainty: The Discounted Case}, + journal = {Journal of Economic Theory}, + volume = {4}, + number = {3}, + pages = {479--513}, + year = {1972} +} + +@article{Sargent1980q, + author = {Sargent, Thomas J.}, + title = {``{T}obin's q'' and the Rate of Investment in General Equilibrium}, + journal = {Carnegie-Rochester Conference Series on Public Policy}, + volume = {12}, + pages = {107--154}, + year = {1980} +} + +@book{Krishna2009, + author = {Krishna, Vijay}, + title = {Auction Theory}, + edition = {2nd}, + publisher = {Academic Press}, + address = {San Diego}, + year = {2009} +} + +@article{ShapleyShubik1971, + author = {Shapley, Lloyd S. and Shubik, Martin}, + title = {The Assignment Game {I}: The Core}, + journal = {International Journal of Game Theory}, + volume = {1}, + number = {1}, + pages = {111--130}, + year = {1971} +} + +@article{CrawfordKnoer1981, + author = {Crawford, Vincent P. and Knoer, Elsie Marie}, + title = {Job Matching with Heterogeneous Firms and Workers}, + journal = {Econometrica}, + volume = {49}, + number = {2}, + pages = {437--450}, + year = {1981} +} + +@article{DemangeGaleSotomayor1986, + author = {Demange, Gabrielle and Gale, David and Sotomayor, Marilda}, + title = {Multi-Item Auctions}, + journal = {Journal of Political Economy}, + volume = {94}, + number = {4}, + pages = {863--872}, + year = {1986} +} + +@article{Leonard1983, + author = {Leonard, Herman B.}, + title = {Elicitation of Honest Preferences for the Assignment of + Individuals to Positions}, + journal = {Journal of Political Economy}, + volume = {91}, + number = {3}, + pages = {461--479}, + year = {1983} +} + +@article{EdelmanOstrovskySchwarz2007, + author = {Edelman, Benjamin and Ostrovsky, Michael and Schwarz, Michael}, + title = {Internet Advertising and the Generalized Second-Price Auction: + Selling Billions of Dollars Worth of Keywords}, + journal = {American Economic Review}, + volume = {97}, + number = {1}, + pages = {242--259}, + year = {2007} +} + +@article{RileySamuelson1981, + author = {Riley, John G. and Samuelson, William F.}, + title = {Optimal Auctions}, + journal = {American Economic Review}, + volume = {71}, + number = {3}, + pages = {381--392}, + year = {1981} +} + +@article{Klemperer1999, + author = {Klemperer, Paul}, + title = {Auction Theory: A Guide to the Literature}, + journal = {Journal of Economic Surveys}, + volume = {13}, + number = {3}, + pages = {227--286}, + year = {1999} +} + +@book{Milgrom2004, + author = {Milgrom, Paul}, + title = {Putting Auction Theory to Work}, + publisher = {Cambridge University Press}, + address = {Cambridge}, + year = {2004} +} + +@book{DavidNagaraja2003, + author = {David, Herbert A. and Nagaraja, Haikady N.}, + title = {Order Statistics}, + edition = {3rd}, + publisher = {Wiley}, + address = {Hoboken, NJ}, + year = {2003} +} + +@techreport{AlvarezArgente2026, + author = {Alvarez, Fernando and Argente, David}, + title = {Strategic Complementarity and Persistence in Linear-Quadratic + Mean Field Games}, + institution = {University of Chicago and Yale University}, + year = {2026} +} + +@article{LasryLions2007, + author = {Lasry, Jean-Michel and Lions, Pierre-Louis}, + title = {Mean Field Games}, + journal = {Japanese Journal of Mathematics}, + volume = {2}, + number = {1}, + pages = {229--260}, + year = {2007} +} + +@article{HuangMalhameCaines2006, + author = {Huang, Minyi and Malham{\'e}, Roland P. and Caines, Peter E.}, + title = {Large Population Stochastic Dynamic Games: Closed-Loop + {McKean-Vlasov} Systems and the Nash Certainty Equivalence + Principle}, + journal = {Communications in Information and Systems}, + volume = {6}, + number = {3}, + pages = {221--252}, + year = {2006} +} + +@article{AchdouEtAl2022, + author = {Achdou, Yves and Han, Jiequn and Lasry, Jean-Michel and + Lions, Pierre-Louis and Moll, Benjamin}, + title = {Income and Wealth Distribution in Macroeconomics: A + Continuous-Time Approach}, + journal = {Review of Economic Studies}, + volume = {89}, + number = {1}, + pages = {45--86}, + year = {2022} +} + +@book{CarmonaDelarue2018, + author = {Carmona, Ren{\'e} and Delarue, Fran{\c c}ois}, + title = {Probabilistic Theory of Mean Field Games with Applications}, + publisher = {Springer}, + address = {Cham}, + year = {2018} +} + +@book{LancasterRodman1995, + author = {Lancaster, Peter and Rodman, Leiba}, + title = {Algebraic {R}iccati Equations}, + publisher = {Oxford University Press}, + address = {Oxford}, + year = {1995} +} + +@article{Kimball1995, + author = {Kimball, Miles S.}, + title = {The Quantitative Analytics of the Basic Neomonetarist Model}, + journal = {Journal of Money, Credit and Banking}, + volume = {27}, + number = {4}, + pages = {1241--1277}, + year = {1995} +} + +@article{KlenowWillis2016, + author = {Klenow, Peter J. and Willis, Jonathan L.}, + title = {Real Rigidities and Nominal Price Changes}, + journal = {Economica}, + volume = {83}, + number = {331}, + pages = {443--472}, + year = {2016} +} + +@article{Midrigan2011, + author = {Midrigan, Virgiliu}, + title = {Menu Costs, Multiproduct Firms, and Aggregate Fluctuations}, + journal = {Econometrica}, + volume = {79}, + number = {4}, + pages = {1139--1180}, + year = {2011} +} + +@article{AlvarezLippi2014, + author = {Alvarez, Fernando and Lippi, Francesco}, + title = {Price Setting with Menu Cost for Multiproduct Firms}, + journal = {Econometrica}, + volume = {82}, + number = {1}, + pages = {89--135}, + year = {2014} +} diff --git a/lectures/_toc.yml b/lectures/_toc.yml index 8fce349a0..953297f68 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -159,6 +159,8 @@ parts: - file: lake_model - file: endogenous_lake - file: rational_expectations + - file: lucas_prescott_investment + - file: optimal_growth_uncertainty - file: re_with_feedback - file: markov_perf - file: uncertainty_traps @@ -167,6 +169,7 @@ parts: - file: ak2 - file: ak_aiyagari - file: two_computation + - file: lq_mean_field_games - caption: Asset Pricing and Finance numbered: true chapters: diff --git a/lectures/hoist_failure.md b/lectures/hoist_failure.md index a991b426c..61290afd5 100644 --- a/lectures/hoist_failure.md +++ b/lectures/hoist_failure.md @@ -276,6 +276,8 @@ print(f"Sample mean: {samp_mean:.3f}") We define helper functions to create discretized versions of lognormal probability density functions. +We write out the density by hand to keep the formula {eq}`lognormal_pdf` in view; `scipy.stats.lognorm(s=σ, scale=np.exp(μ)).pdf(x)` computes the same thing. + ```{code-cell} ipython3 def lognormal_pdf(x, μ, σ): """ @@ -288,7 +290,19 @@ def lognormal_pdf(x, μ, σ): def discretize_lognormal(μ, σ, I, m): """ - Create discretized lognormal probability mass function. + Discretize a lognormal distribution on the grid 0, m, 2m, ..., up to I. + + Parameters + ---------- + μ, σ : parameters of the lognormal distribution + I : upper end of the grid, which truncates the right tail + m : spacing between grid points, which sets the resolution + + Returns + ------- + p_array : the density evaluated on the grid + p_array_norm : the implied probability mass function, summing to one + x : the grid itself, with I / m points """ x = np.arange(1e-7, I, m) p_array = lognormal_pdf(x, μ, σ) @@ -296,17 +310,26 @@ def discretize_lognormal(μ, σ, I, m): return p_array, p_array_norm, x ``` -We set the grid length $I$ to a power of 2 to enable efficient Fast Fourier Transform computation. +Two separate choices govern the quality of this approximation, and it pays to keep them straight. + +* $I$ fixes where the grid **stops**, so it controls how much of the right tail we throw away +* $m$ fixes the **spacing** between grid points, so it controls resolution + +The grid has $I/m$ points, so raising $I$ at fixed $m$ buys range, while lowering $m$ at fixed $I$ buys accuracy. + +Once $I$ is large enough that almost no probability mass lies beyond it, further increases change nothing, and only $m$ matters. + +{ref}`hoist_ex1` asks you to verify this. ```{note} -Increasing the power $p$ (e.g., from 12 to 15) improves the approximation quality but increases computational cost. +`scipy.signal.fftconvolve` pads its inputs to a convenient length internally, so there is no need to choose $I/m$ to be a power of two. ``` ```{code-cell} ipython3 # Set grid parameters p = 15 -I = 2**p # Truncation value (power of 2 for FFT efficiency) -m = 0.1 # Increment size +I = 2**p # where the grid stops: truncates the right tail +m = 0.1 # spacing between grid points: sets the resolution ``` Let's visualize how well the discretized distribution approximates the continuous lognormal distribution: @@ -347,7 +370,7 @@ print(f"Discretized mean: {mean_discrete:.3f}") Now let's use the convolution theorem to compute the probability distribution of a sum of the two lognormal random variables we have parameterized above. -We'll also compute the probability of a sum of three log normal distributions constructed above. +We'll also compute the probability distribution of a sum of three log normal distributions constructed above. For long sequences, `scipy.signal.fftconvolve` is much faster than `numpy.convolve` because it uses Fast Fourier Transforms. @@ -390,30 +413,39 @@ This is the algorithm used by `fftconvolve`. Let's do a warmup calculation that compares the times taken by `numpy.convolve` and `scipy.signal.fftconvolve` +Our three components are identically distributed, so a single discretization serves for all of them. + ```{code-cell} ipython3 -# Discretize three lognormal distributions +# Discretize the lognormal distribution; the three components are i.i.d. _, pmf1, x = discretize_lognormal(μ, σ, I, m) -_, pmf2, x = discretize_lognormal(μ, σ, I, m) -_, pmf3, x = discretize_lognormal(μ, σ, I, m) +pmf2 = pmf3 = pmf1 + +# Direct convolution costs O(N²), so we time it on a short prefix +short = pmf1[:20_000] -# Time numpy.convolve with qe.Timer() as timer_numpy: - conv_np = np.convolve(pmf1, pmf2) - conv_np = np.convolve(conv_np, pmf3) + np.convolve(short, short) time_numpy = timer_numpy.elapsed -# Time fftconvolve with qe.Timer() as timer_fft: - conv_fft = fftconvolve(pmf1, pmf2) - conv_fft = fftconvolve(conv_fft, pmf3) + fftconvolve(short, short) time_fft = timer_fft.elapsed -print(f"Time with np.convolve: {time_numpy:.4f} seconds") -print(f"Time with fftconvolve: {time_fft:.4f} seconds") -print(f"Speedup factor: {time_numpy / time_fft:.1f}x") +print(f"On {len(short):,} points:") +print(f" np.convolve: {time_numpy:.4f} seconds") +print(f" fftconvolve: {time_fft:.4f} seconds") +print(f" speedup: {time_numpy / time_fft:.0f}x") ``` -The Fast Fourier Transform provides orders of magnitude speedup. +The gap widens rapidly with the length of the sequences, because direct convolution costs $O(N^2)$ operations while the FFT approach costs $O(N \log N)$. + +On the full grid used below, the direct method is slower by more than two orders of magnitude. + +```{code-cell} ipython3 +# The full calculation, done the fast way +conv_fft = fftconvolve(fftconvolve(pmf1, pmf2), pmf3) +print(f"grid points per component: {len(pmf1):,}") +``` Now let’s plot our computed probability mass function approximation for the sum of two log normal random variables against the histogram of the sample that we formed above @@ -482,7 +514,7 @@ print(f" Computed mean: {mean_conv3:.3f}") We shall soon apply the convolution theorem to compute the probability of a **top event** in a failure tree analysis. -Before applying the convolution theorem, we first describe the model that connects constituent events to the *top* end whose failure rate we seek to quantify. +Before applying the convolution theorem, we first describe the model that connects constituent events to the *top event* whose failure rate we seek to quantify. Fault tree analysis is a widely used technique for assessing system reliability, as described by {cite:t}`Ardron_2018`. @@ -521,7 +553,7 @@ We assume: * The failure probability $P(A_i)$ of each component $A_i$ is small * Component failures are statistically independent -We repeatedly apply a **rare event approximation** to obtain the following formula for the problem of a system failure: +We repeatedly apply a **rare event approximation** to obtain the following formula for the probability of a system failure: $$ P(F) \approx P(A_1) + P (A_2) + \cdots + P (A_n) @@ -539,6 +571,14 @@ where $P(F)$ is the system failure probability. Probabilities for each event are recorded as failure rates per year. +```{note} +Strictly speaking, a failure *rate* per year and a failure *probability* within a year are different objects. + +For rare events they nearly coincide, because $1 - e^{-\lambda} \approx \lambda$ when $\lambda$ is small. + +The same approximation that lets us add probabilities across components also lets us move between rates and probabilities, so we follow the reliability literature in using the two words interchangeably here. +``` + ## Failure rates unknown Now we come to the problem that really interests us, following {cite:t}`Ardron_2018` and @@ -551,7 +591,7 @@ We address this problem by specifying **probabilities of probabilities** that c Thus, we assume that a system analyst is uncertain about the failure rates $P(A_i), i =1, \ldots, n$ for components of a system. -The analyst copes with this situation by regarding the systems failure probability $P(F)$ and each of the component probabilities $P(A_i)$ as random variables. +The analyst copes with this situation by regarding the system's failure probability $P(F)$ and each of the component probabilities $P(A_i)$ as random variables. * dispersions of the probability distribution of $P(A_i)$ characterizes the analyst's uncertainty about the failure probability $P(A_i)$ @@ -559,6 +599,18 @@ The analyst copes with this situation by regarding the systems failure probabili This leads to what is sometimes called a **hierarchical** model in which the analyst has probabilities about the probabilities $P(A_i)$. +```{note} +Two distinct kinds of randomness appear in this model, and it is worth keeping them apart. + +*Aleatory* uncertainty is the randomness in whether a component fails during a given year; it is described by the failure rate $P(A_i)$. + +*Epistemic* uncertainty is the analyst's ignorance about the value of that rate; it is described by the lognormal distribution that he places over $P(A_i)$. + +The distribution that we compute below is an epistemic object: it describes what the analyst knows about a failure rate, not how often the system fails. + +Separating the two is the central recommendation of {cite:t}`apostolakis1990`. +``` + The analyst formalizes his uncertainty by assuming that * the failure probability $P(A_i)$ is itself a log normal random variable with parameters $(\mu_i, \sigma_i)$. @@ -570,8 +622,18 @@ The analyst assumes that such information about the observed dispersion of annu The analyst assumes that the random variables $P(A_i)$ are statistically mutually independent. +```{warning} +Independence is a strong assumption and it is the one that reliability analysts worry about most. + +A design flaw, a shared power supply, a common maintenance crew, or a single environmental shock can push many components toward failure at once. + +Such **common-cause** failures make the upper tail of the distribution of $P(F)$ much fatter than the independent calculation suggests, which is precisely the region that a safety regulator cares about. + +{ref}`hoist_ex5` quantifies how much difference this makes. +``` + The analyst wants to approximate a probability mass function and cumulative distribution function -of the systems failure probability $P(F)$. +of the system's failure probability $P(F)$. * We say probability mass function because of how we discretize each random variable, as described earlier. @@ -588,12 +650,6 @@ A regulatory agency requires the system to be designed so that the top event fai ### Model specification -We'll take close to a real world example by assuming that $n = 14$. - -The example estimates the annual failure rate of a critical hoist at a nuclear waste facility. - -A regulatory agency wants the system to be designed in a way that makes the failure rate of the top event small with high probability. - This example is Design Option B-2 (Case I) described in Table 10 on page 27 of {cite:t}`Greenfield_Sargent_1993`. The table describes parameters $\mu_i, \sigma_i$ for fourteen log normal random variables that consist of **seven pairs** of random variables that are identically and independently distributed. @@ -628,7 +684,12 @@ We define a helper function to find array indices: ```{code-cell} ipython3 def find_nearest(array, value): """ - Find the index of the array element nearest to the given value. + Index of the array element nearest to the given value. + + Applied to a cumulative distribution function, this returns the grid point + whose cumulative probability is closest to a target, which for a finely + discretized distribution is indistinguishable from the usual definition of + a quantile as the smallest x with CDF(x) >= q. """ array = np.asarray(array) idx = (np.abs(array - value)).argmin() @@ -665,8 +726,33 @@ with qe.Timer() as timer: system_pmf = fftconvolve(system_pmf, pmf) print(f"Time for 13 convolutions: {timer.elapsed:.4f} seconds") + +# the convolution lives on the same grid spacing, but extends much further +system_grid = np.arange(len(system_pmf)) * m +print(f"grid points in the answer: {len(system_pmf):,}") +``` + +Before plotting the cumulative distribution function, let's look at the density itself. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Density of the system failure rate + name: fig-hoist-pdf +--- +fig, ax = plt.subplots(figsize=(10, 6)) +upper = 2000 +ax.plot(system_grid[:int(upper/m)], system_pmf[:int(upper/m)] / m, 'b-', lw=2) +ax.set_xlabel(r'failure rate ($\times 10^{-9}$ per year)') +ax.set_ylabel('density') +plt.show() ``` +The density is strongly skewed to the right: a long upper tail stretches far beyond the bulk of the distribution. + +This asymmetry is what makes a single point estimate of a failure rate a poor summary, and it is why the analyst reports quantiles instead. + We now plot a counterpart to the cumulative distribution function (CDF) in figure 5 on page 29 of {cite:t}`Greenfield_Sargent_1993` ```{code-cell} ipython3 @@ -700,136 +786,373 @@ We also present a counterpart to their Table 11 on page 28 of {cite:t}`Greenfiel ```{code-cell} ipython3 -# Find quantiles -quantiles = [0.01, 0.05, 0.10, 0.50, 0.665, 0.85, 0.90, 0.95, 0.99, 0.9978] -quantile_values = [x[find_nearest(cdf, q)] for q in quantiles] +# Percentiles reported in Table 11 of Greenfield and Sargent (1993), +# together with their published values, in units of 10^-9 per year +reference = {1.0: 77, 10.0: 130, 50.0: 263, 66.5: 341, + 85.0: 513, 95.0: 811, 99.0: 1480, 99.78: 2490} -# Create table -table_data = [[f"{100*q:.2f}%", f"{val:.3f}"] - for q, val in zip(quantiles, quantile_values)] +table_data = [] +for pc, published in reference.items(): + ours = system_grid[find_nearest(cdf, pc/100)] + table_data.append([f"{pc}%", f"{ours:.1f}", published, + f"{100*(ours - published)/published:+.1f}%"]) print("\nSystem failure rate quantiles (×10^-9 per year):") -print(tabulate(table_data, - headers=['Percentile', 'Failure rate'], tablefmt='grid')) +print(tabulate(table_data, + headers=['Percentile', 'Computed here', 'Greenfield-Sargent', 'Difference'], + tablefmt='grid')) +``` + +Our quantiles reproduce the published ones to within one and a half per cent, and slightly understate each of them. + +The small discrepancies reflect the precision of the reported parameters $\mu_i, \sigma_i$, the grid spacing $m$, and the point at which the grid is truncated. + +### Reading the answer + +The numbers in this table, rather than any single one of them, are the output of the analysis. + +The median failure rate is about $261 \times 10^{-9}$ per year, while the 95th percentile is about $808 \times 10^{-9}$, three times larger. + +That spread is not a statement about how often the hoist fails; it is a statement about how little the analyst knows about how often the hoist fails. + +Notice also where the *mean* of the distribution falls. + +```{code-cell} ipython3 +mean_rate = np.sum(system_grid * system_pmf) +mean_percentile = 100 * cdf[find_nearest(system_grid, mean_rate)] + +print(f"mean failure rate: {mean_rate:.1f} × 10⁻⁹ per year") +print(f"the mean sits at the {mean_percentile:.1f}th percentile") +``` + +Because the distribution is skewed, the mean lies well above the median, at about the 66th percentile. + +This is why Table 11 of {cite:t}`Greenfield_Sargent_1993` records the mean alongside the 66.5th percentile. + +The practical implication is the one that motivated the original study. + +```{code-cell} ipython3 +# The point estimate used in the U.S. Department of Energy's 1990 risk +# assessment, expressed in the same units +doe_estimate = 220 # 2.2 × 10^-7 per year + +pct = 100 * cdf[find_nearest(system_grid, doe_estimate)] +print(f"the DOE point estimate of {doe_estimate} × 10⁻⁹ lies at " + f"the {pct:.0f}th percentile") +print(f"so the analyst assigns probability {100-pct:.0f}% to the " + f"true rate exceeding it") ``` -The computed quantiles agree closely with column 2 of Table 11 on page 28 of {cite}`Greenfield_Sargent_1993`. +An analysis that reports a single number in place of a distribution conveys none of this. -Minor discrepancies may be due to differences in: -* Numerical precision of input parameters $\mu_i, \sigma_i$ -* Number of grid points in the discretization -* Grid increment size +{cite:t}`Greenfield_Sargent_1993` made exactly this point: reading their figure, they put the Department of Energy's point estimate at the 36th percentile and concluded that there was roughly a 64 per cent chance that the true failure rate was higher. ## Exercises -```{exercise-start} +```{exercise} :label: hoist_ex1 -``` -Experiment with different values of the power parameter $p$ (which determines the grid size as $I = 2^p$). +Our discretization involves two separate choices: where the grid stops, $I = 2^p$, and how finely it is spaced, $m$. -Try $p \in \{12, 13, 14, 15, 16\}$ and compare: -1. Computation time -2. Accuracy of the median (50th percentile) compared to the reference value -3. Memory usage implications +Investigate what each one controls. -What trade-offs do you observe? -```{exercise-end} +1. Holding $m = 0.05$ fixed, compute the median, the 95th percentile and the 99.78th percentile of the system failure rate for $p = 10, 11, \ldots, 15$. For each $p$, also compute how much probability mass the truncation discards, using $\sum_i \Pr(P(A_i) > I)$. +1. Holding $p = 14$ fixed, repeat for $m = 0.4, 0.2, 0.1, 0.05, 0.025$. +1. Which statistic is sensitive to which choice, and why? Are the values $p = 15$, $m = 0.05$ used in the lecture well chosen? ``` ```{solution-start} hoist_ex1 :class: dropdown ``` -Here is one solution: - ```{code-cell} ipython3 -# Test different grid sizes -p_values = [12, 13, 14, 15, 16] -results = [] +from scipy.stats import norm + +def system_distribution(p_grid, m_grid): + "Failure rate distribution of the whole system on a given grid." + I_grid = 2**p_grid + pmfs = [] + for μ_i, σ_i in params[:6]: + _, pmf_i, _ = discretize_lognormal(μ_i, σ_i, I_grid, m_grid) + pmfs.append(pmf_i) + μ7, σ7 = params[6] + _, pmf7, _ = discretize_lognormal(μ7, σ7, I_grid, m_grid) + pmfs.extend([pmf7] * 8) + + total = pmfs[0] + for pmf_i in pmfs[1:]: + total = fftconvolve(total, pmf_i) + return total, np.arange(len(total)) * m_grid + -for p_test in p_values: - I_test = 2**p_test - m_test = 0.05 +def quantiles_of(pmf, grid, levels=(0.5, 0.95, 0.9978)): + cdf_local = np.cumsum(pmf) + return [grid[find_nearest(cdf_local, q)] for q in levels] - # Discretize distributions - pmfs_test = [] - for μ, σ in params[:6]: - _, pmf, x_test = discretize_lognormal(μ, σ, I_test, m_test) - pmfs_test.append(pmf) - # Add 8 copies of component type 7 +def discarded_mass(I_grid): + "Probability that a component's rate exceeds the end of the grid." + lost = sum(norm.sf((np.log(I_grid) - μ_i)/σ_i) for μ_i, σ_i in params[:6]) μ7, σ7 = params[6] - _, pmf7, x_test = discretize_lognormal(μ7, σ7, I_test, m_test) - pmfs_test.extend([pmf7] * 8) + return lost + 8 * norm.sf((np.log(I_grid) - μ7)/σ7) - # Time the convolutions - with qe.Timer() as timer_test: - system_test = pmfs_test[0] - for pmf in pmfs_test[1:]: - system_test = fftconvolve(system_test, pmf) - # Compute median - cdf_test = np.cumsum(system_test) - median = x_test[find_nearest(cdf_test, 0.5)] +rows = [] +for p_test in range(10, 16): + pmf_t, grid_t = system_distribution(p_test, 0.05) + med, q95, q9978 = quantiles_of(pmf_t, grid_t) + rows.append([p_test, 2**p_test, f"{discarded_mass(2**p_test):.1e}", + f"{med:.2f}", f"{q95:.2f}", f"{q9978:.2f}"]) - results.append([p_test, I_test, - f"{timer_test.elapsed:.4f}", f"{median:.7f}"]) +print(tabulate(rows, headers=['p', 'I', 'mass discarded', + 'median', '95th', '99.78th'], tablefmt='grid')) +``` -print(tabulate(results, - headers=['p', 'Grid size (2^p)', 'Time (s)', 'Median'], +```{code-cell} ipython3 +rows = [] +for m_test in (0.4, 0.2, 0.1, 0.05, 0.025): + pmf_t, grid_t = system_distribution(14, m_test) + med, q95, q9978 = quantiles_of(pmf_t, grid_t) + rows.append([m_test, len(grid_t), f"{med:.3f}", f"{q95:.2f}", f"{q9978:.2f}"]) + +print(tabulate(rows, headers=['m', 'grid points', 'median', '95th', '99.78th'], tablefmt='grid')) ``` -The results typically show the following trade-offs: -- Larger grid sizes provide better accuracy but increase computation time -- The relationship between $p$ and computation time is roughly linear for FFT-based convolution -- Beyond $p = 13$, the accuracy gains diminish while computational cost continues to grow -- For this application, $p = 13$ provides a good balance between accuracy and efficiency +The two choices do quite different jobs. + +Truncation governs the **far tail**. At $p = 10$ the 99.78th percentile is badly understated, and it keeps rising until about $p = 14$, by which point the discarded mass has fallen to roughly $10^{-6}$. The median, by contrast, has settled by $p = 12$: throwing away the extreme right tail of each component hardly moves the middle of the distribution of their sum. + +Resolution governs **overall precision**. Halving $m$ shifts every quantile slightly and uniformly, and the shifts are small: going from $m = 0.4$ to $m = 0.025$ moves the median by about 1.5 per cent. + +The lecture's choices are sensible. With $p = 15$ the discarded mass is around $10^{-7}$, so even the 99.78th percentile is accurate, and $m = 0.05$ is fine enough that further refinement changes little. + +The moral is that a grid that looks adequate for the median can be badly inadequate for the upper tail, which is exactly the region a safety regulator cares about. ```{solution-end} ``` -```{exercise-start} +```{exercise} :label: hoist_ex2 + +The rare event approximation replaces $P(A \cup B)$ by $P(A) + P(B)$, discarding $P(A \cap B)$. + +Assess how good it is here. + +1. Using the *mean* failure rate of each of the fourteen components as a representative value, compare $\sum_i p_i$ with the exact probability that at least one component fails, $1 - \prod_i (1 - p_i)$. +1. Repeat with all the rates multiplied by $10^3$, $10^6$ and $10^7$, and report the relative error in each case. +1. At what order of magnitude does the approximation start to matter? ``` -The rare event approximation assumes that $P(A_i) P(A_j)$ is negligible compared to $P(A_i) + P(A_j)$. +```{solution-start} hoist_ex2 +:class: dropdown +``` -Using the computed distribution, calculate the expected value of the system failure rate and compare it to the sum of the expected values of the individual component failure rates. +```{code-cell} ipython3 +# representative rate for each of the 14 components +component_means = [np.exp(μ_i + 0.5*σ_i**2) for μ_i, σ_i in params[:6]] +μ7, σ7 = params[6] +component_means.extend([np.exp(μ7 + 0.5*σ7**2)] * 8) +component_means = np.array(component_means) + +rows = [] +for factor, label in ((1e-9, 'as calibrated'), (1e-6, '× 10³'), + (1e-3, '× 10⁶'), (1e-2, '× 10⁷')): + probs = component_means * factor + approx = probs.sum() + exact = 1 - np.prod(1 - probs) + rows.append([label, f"{approx:.6e}", f"{exact:.6e}", + f"{100*(approx - exact)/exact:.4f}%"]) + +print(tabulate(rows, headers=['failure rates', 'Σ pᵢ', '1 - Π(1-pᵢ)', + 'relative error'], tablefmt='grid')) +``` -How good is the rare event approximation in this case? -```{exercise-end} +At the calibrated magnitudes, around $3 \times 10^{-7}$ per year in total, the approximation is exact to the precision shown: the neglected term is of order $p_i p_j \approx 10^{-14}$. + +Multiplying every rate by a thousand still leaves an error of about one part in ten thousand. + +The approximation only becomes consequential once individual failure probabilities reach the order of a per cent, where it overstates the system failure probability by more than ten per cent, and it fails completely when $\sum_i p_i$ approaches or exceeds one, where it can return a "probability" greater than one. + +Note that a naive check of this approximation -- comparing the mean of the computed distribution of $\sum_i P(A_i)$ with the sum of the component means -- reveals nothing, because those two quantities are equal by linearity of expectation whatever the quality of the approximation. + +```{solution-end} ``` +```{exercise} +:label: hoist_ex3 -```{solution-start} hoist_ex2 +A regulator who learns that the 95th percentile of the failure rate is too high will ask which components to improve. + +Answer that question by computing, for each of the seven component types, the 95th percentile of the system failure rate when that type is removed from the system entirely. + +Rank the component types by how much they contribute to the upper tail, and compare the ranking with the components' mean failure rates. +``` + +```{solution-start} hoist_ex3 :class: dropdown ``` -Here is one solution: +```{code-cell} ipython3 +def system_without(drop): + "System failure rate distribution with component type `drop` removed." + pmfs = [] + for k, (μ_i, σ_i) in enumerate(params[:6]): + if k == drop: + continue + _, pmf_i, _ = discretize_lognormal(μ_i, σ_i, I, m) + pmfs.append(pmf_i) + if drop != 6: + μ7, σ7 = params[6] + _, pmf7, _ = discretize_lognormal(μ7, σ7, I, m) + pmfs.extend([pmf7] * 8) + + total = pmfs[0] + for pmf_i in pmfs[1:]: + total = fftconvolve(total, pmf_i) + return total, np.arange(len(total)) * m + + +base_q95 = system_grid[find_nearest(cdf, 0.95)] + +rows = [] +for k in range(7): + pmf_k, grid_k = system_without(k) + q95 = grid_k[find_nearest(np.cumsum(pmf_k), 0.95)] + μ_k, σ_k = params[k] + n_units = 8 if k == 6 else 1 + rows.append([f"type {k+1}", n_units, f"{np.exp(μ_k + 0.5*σ_k**2):.1f}", + f"{q95:.1f}", f"{100*(base_q95 - q95)/base_q95:.1f}%"]) + +rows.sort(key=lambda r: -float(r[4].rstrip('%'))) +print(f"95th percentile with all components: {base_q95:.1f}\n") +print(tabulate(rows, headers=['removed', 'units', 'mean rate each', + '95th pct without it', 'reduction'], + tablefmt='grid')) +``` + +Component type 1 dominates: removing that single unit cuts the 95th percentile by almost half, far more than any other change available to the designer. + +The ranking follows the components' mean rates closely here, because all seven types have similar dispersions. + +It need not do so in general: a component with a modest mean but a large $\sigma$ contributes disproportionately to the upper tail, which is why the analyst works with the whole distribution rather than with means. + +Note also that type 7, which appears eight times, matters less than type 1, which appears once. + +Counting components is no guide to where the risk lies. + +```{solution-end} +``` + +```{exercise} +:label: hoist_ex4 + +We could have computed the distribution of the system failure rate by simulation instead of by convolution. + +Draw samples of all fourteen component rates, sum them, and compare the resulting quantiles with those from the convolution, for sample sizes $10^4$, $10^5$ and $10^6$. + +Compare the median, the 95th percentile and the 99.78th percentile. + +Which method would you prefer, and why? +``` + +```{solution-start} hoist_ex4 +:class: dropdown +``` ```{code-cell} ipython3 -# Create extended grid for the convolution result -x_extended = np.arange(0, len(system_pmf) * m, m) -E_system = np.sum(x_extended * system_pmf) +all_params = list(params[:6]) + [params[6]] * 8 +rng_mc = np.random.default_rng(0) +levels = (50, 95, 99.78) -# Compute sum of individual expected values -component_means = [np.exp(μ + 0.5 * σ**2) for μ, σ in params[:6]] -# Add 8 components of type 7 -μ7, σ7 = params[6] -component_means.extend([np.exp(μ7 + 0.5 * σ7**2)] * 8) +rows = [] +for N in (10_000, 100_000, 1_000_000): + draws = sum(rng_mc.lognormal(μ_i, σ_i, N) for μ_i, σ_i in all_params) + rows.append([f"{N:,}"] + [f"{np.percentile(draws, pc):.1f}" for pc in levels]) + +rows.append(['convolution'] + + [f"{system_grid[find_nearest(cdf, pc/100)]:.1f}" for pc in levels]) + +print(tabulate(rows, headers=['method', 'median', '95th', '99.78th'], + tablefmt='grid')) +``` + +Simulation converges to the same answer, which is a useful check on both calculations. + +The two methods differ in where their errors lie. + +Monte Carlo error is largest exactly where the analysis matters most: the 99.78th percentile is pinned down by roughly one draw in five hundred, so with $10^4$ draws only about twenty observations inform it, and the estimate is visibly off. + +The convolution, by contrast, computes the whole distribution at once and its error comes from the grid rather than from sampling noise, so it is equally accurate in the tail as in the middle. + +It is also far quicker: thirteen fast convolutions take a couple of seconds, and the answer does not change when you rerun it with a different seed. + +```{solution-end} +``` + +```{exercise} +:label: hoist_ex5 + +The entire calculation assumes that the fourteen component failure rates are statistically independent. + +Investigate what happens when they are not. + +Suppose that -E_sum = sum(component_means) +$$ +\log P(A_i) = \mu_i + \sigma_i \left( \sqrt{\rho}\, z_0 + \sqrt{1-\rho}\, z_i \right), +$$ + +where $z_0$ is a shock common to all components and $z_1, \ldots, z_{14}$ are idiosyncratic, all standard normal. + +Each component still has exactly its original marginal distribution, but any two of them now have correlation $\rho$ in logs. + +Simulate the system failure rate for $\rho = 0, 0.2, 0.5, 0.8$ and report the median, the 95th, the 99th and the 99.9th percentiles. -print(f"Expected system failure rate: {E_system:.3f} × 10^-9") -print(f"Sum of component expected failure rates: {E_sum:.3f} × 10^-9") -print(f"Relative difference: {100 * abs(E_system - E_sum) / E_sum:.2f}%") +Explain what happens and why it matters for a safety analysis. ``` -The rare event approximation works well when failure probabilities are small. +```{solution-start} hoist_ex5 +:class: dropdown +``` + +```{code-cell} ipython3 +μ_vec = np.array([q[0] for q in all_params]) +σ_vec = np.array([q[1] for q in all_params]) + +N_sim = 400_000 +rng_cc = np.random.default_rng(1) + +rows = [] +for ρ in (0.0, 0.2, 0.5, 0.8): + z0 = rng_cc.normal(size=(N_sim, 1)) + zi = rng_cc.normal(size=(N_sim, len(all_params))) + logs = μ_vec + σ_vec * (np.sqrt(ρ)*z0 + np.sqrt(1-ρ)*zi) + totals = np.exp(logs).sum(axis=1) + rows.append([ρ] + [f"{np.percentile(totals, pc):.0f}" + for pc in (50, 95, 99, 99.9)]) + +print(tabulate(rows, headers=['ρ', 'median', '95th', '99th', '99.9th'], + tablefmt='grid')) +``` + +Correlation leaves each component's marginal distribution untouched, and it leaves the mean of the sum untouched as well. + +What it changes is the shape of the distribution of the sum. + +With independent components, a high draw for one is typically offset by ordinary draws for the others, and the fourteen-fold averaging produces a relatively concentrated total. + +A common shock removes that diversification: when $z_0$ is large every component is bad at once. + +The result is a distribution with a *lower* median and a much *heavier* upper tail. + +At $\rho = 0.8$ the median falls by about a third while the 99.9th percentile rises by roughly three quarters. + +For a safety analysis this is the dangerous direction of error. + +Assuming independence when a common cause is present makes the system look both typically safer and much less likely to suffer a very bad year than it really is. -The expected value of the sum equals the sum of the expected values (by linearity of expectation), so these should match closely regardless of the rare event approximation. +This is why reliability studies devote so much attention to identifying shared power supplies, shared maintenance procedures, common design faults, and other mechanisms that defeat the independence assumption. ```{solution-end} ``` diff --git a/lectures/house_auction.md b/lectures/house_auction.md index 6dc8eb9a7..5f8a1cca3 100644 --- a/lectures/house_auction.md +++ b/lectures/house_auction.md @@ -41,7 +41,7 @@ We describe two distinct mechanisms * A multiple rounds, ascending bid auction - * A special case of a Groves-Clarke {cite}`Groves_73`, {cite}`Clarke_71` mechanism with a benevolent social planner + * A sequential, "greedy" second-price mechanism administered by a benevolent social planner, inspired by the ideas of {cite}`Groves_73` and {cite}`Clarke_71` ```{note} In 1994, the multiple rounds, ascending bid auction was actually used by Stanford University to sell leases to 9 lots on the Stanford campus to eligible faculty members. @@ -77,7 +77,7 @@ There are multiple **rounds** - all $n$ houses are allocated after the final round - - house $i$ is retained by the auctioneer if not prospective buyer offers more that $r_i$ for the house + - house $i$ is retained by the auctioneer if no prospective buyer offers more than $r_i$ for the house In this auction, person $j$ never tells anyone else his/her private values $v_{ij}$ @@ -92,18 +92,25 @@ Among all feasible allocations, a **socially optimal allocation** maximizes the The planner tells everyone in advance how he/she will allocate houses based on the matrix of values that prospective buyers report. -The mechanism provide every prospective buyer an incentive to reveal his vector of private values to the planner. +The mechanism is intended to give every prospective buyer an incentive to reveal his vector of private values to the planner. After the planner receives everyone's vector of private values, the planner deploys a **sequential** algorithm to determine an **allocation** of houses and a set of **fees** that he charges awardees for the negative **externality** that their presence impose on other prospective buyers. -## Equivalence of Allocations +## Comparing the Two Mechanisms Remarkably, these two mechanisms can produce virtually identical allocations. -We construct Python code for both mechanism. +But "virtually identical" is not "identical", and the differences turn out to be instructive. + +We'll discover that + +* the ascending bid auction maximizes the total value of the allocation, and, as the bid increment $\epsilon$ becomes small, the prices that it sets converge to the payments prescribed by the **Vickrey-Clarke-Groves** (VCG) mechanism +* the greedy mechanism of the benevolent planner does **not** always maximize total value, and the prices that it charges are **not** VCG payments, so truthful reporting is not always in a buyer's interest + +We construct Python code for both mechanisms. We also work out some examples by hand or almost by hand. @@ -201,24 +208,34 @@ person $j$ bought house $i$, in which case $S_{ij} = v_{ij} - p_i$ **Proposed buyer strategy:** -In this pseudo code and the actual Python code below, we'll assume that all buyers choose to use the following strategy - - * The strategy is optimal for each buyer +In this pseudo code and the actual Python code below, we'll assume that all buyers choose to use the following strategy, often called **straightforward bidding**: in each round, bid on the house that currently offers the buyer the largest surplus. Each buyer $j = 1, \ldots, m$ uses the same strategy. +Straightforward bidding is *myopic*: a buyer looks only at prevailing prices and ignores how his bid today might affect prices tomorrow. + +That myopia is less restrictive than it might sound. + +Our setting is an **assignment market** of the kind studied by {cite}`ShapleyShubik1971`: each buyer wants at most one house, and values are private. + +For such markets, {cite}`CrawfordKnoer1981` and {cite}`DemangeGaleSotomayor1986` showed that an ascending auction like ours, in which all buyers bid straightforwardly, converges to a competitive equilibrium of the assignment market. + +It converges in particular to the **minimum** competitive equilibrium price vector, the one most favorable to buyers. + +We'll return below to a striking consequence of this fact. + The strategy has the form: - Let $\check p^t$ be the $n \times 1$ vector of prevailing highest-bid prices at the beginning of round $t$ - Let $\epsilon>0$ be the minimum bid increment specified by the seller - For each prospective buyer $j$, compute the index of the best house to bid on during round $t$, namely $\hat i_t = \textrm{argmax}_i\{ [ v_{ij} - \check p^t_i - \epsilon ]\}$ -- If $\max_i\{ [ v_{ij} - \check p^t_i - \epsilon ]\} $ $\leq$ $0$, person $j$ permanently drops out of the auction at round $t$ -- If $v_{\hat i_t, j} - \check p^t_i - \epsilon>0$, person $j$ bids $\check p^t_i + \epsilon$ on house $j$ +- If $\max_i\{ [ v_{ij} - \check p^t_i - \epsilon ]\} \leq 0$, person $j$ permanently drops out of the auction at round $t$ +- If $v_{\hat i_t, j} - \check p^t_{\hat i_t} - \epsilon>0$, person $j$ bids $\check p^t_{\hat i_t} + \epsilon$ on house $\hat i_t$ **Resolving ambiguities**: The protocols we have described so far leave open two possible sources of ambiguity. -(1) **The optimal bid choice for buyers in each round.** It is possible that a buyer has the same surplus value for multiple houses. The argmax function in Python always returns the first argmax element. We instead prefer to randomize among such winner. For that reason, we write our own argmax function below. +(1) **The optimal bid choice for buyers in each round.** It is possible that a buyer has the same surplus value for multiple houses. The argmax function in Python always returns the first argmax element. We instead prefer to randomize among such houses. For that reason, we write our own argmax function below. (2) **Seller's choice of winner if same price bid cast by several buyers.** To resolve this ambiguity, we use the np.random.choice function below. @@ -281,8 +298,8 @@ r ```{code-cell} ipython3 def find_argmax_with_randomness(v): """ - We build our own verion of argmax function such that the argmax index will be returned randomly - when there are multiple maximum values. This function is similiar to np.argmax(v,axis=0) + We build our own version of argmax function such that the argmax index will be returned randomly + when there are multiple maximum values. This function is similar to np.argmax(v,axis=0) Parameters: ---------- @@ -332,10 +349,12 @@ def check_kick_off_condition(v, r, ϵ): r: the reservation price - ϵ: the minimun price increment in each round + ϵ: the minimum price increment in each round """ + n, m = v.shape + # we convert the price vector to a matrix in the same shape as value matrix to facilitate subtraction p_start = (ϵ+r)[:,None] @ np.ones(m)[None,:] @@ -365,7 +384,7 @@ def submit_initial_bid(p_initial, ϵ, v): v: the value matrix - ϵ: the minimun price increment in each round + ϵ: the minimum price increment in each round Returns: ---------- @@ -375,11 +394,14 @@ def submit_initial_bid(p_initial, ϵ, v): """ + n, m = v.shape + buyer_list = np.arange(m) + p = p_initial.copy() p_start_mat = (ϵ + p)[:,None] @ np.ones(m)[None,:] surplus_value = v - p_start_mat - # we only care about active buyers who have positve surplus values + # we only care about active buyers who have positive surplus values active_buyer_diagnosis = (surplus_value > 0).any(axis = 0) active_buyer_list = buyer_list[active_buyer_diagnosis] active_buyer_surplus_value = surplus_value[:,active_buyer_diagnosis] @@ -428,7 +450,7 @@ In this case, buyer 0 has to bid one more time with a higher price, which crowds Therefore, final price could be 3 or 4, depending on the winner in the last round. ```{code-cell} ipython3 -def check_terminal_condition(bid_info, p, v): +def check_terminal_condition(bid_info, p, v, ϵ): """ A function that checks whether the auction ends. @@ -443,9 +465,11 @@ def check_terminal_condition(bid_info, p, v): v: value matrix + ϵ: the minimum price increment in each round + Returns: ---------- - allocation: a dictionary that descirbe how the houses bid are assigned. + allocation: a dictionary that describes how the houses bid are assigned. winner_list: a list of winners @@ -453,6 +477,9 @@ def check_terminal_condition(bid_info, p, v): """ + n, m = v.shape + buyer_list = np.arange(m) + # there may be several buyers bidding one house, we choose a winner randomly winner_list=[np.random.choice(bid_info[ii]) for ii in bid_info.keys()] @@ -470,13 +497,13 @@ def check_terminal_condition(bid_info, p, v): loser_surplus_value = v[:,loser_list] - p_mat loser_decision = (loser_surplus_value > 0).any(axis = 0) - print(~(loser_decision.any())) + print('auction ends:', not loser_decision.any()) return allocation,winner_list,loser_list ``` ```{code-cell} ipython3 -allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v) +allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v, ϵ) ``` ```{code-cell} ipython3 @@ -567,81 +594,35 @@ present_dict(bid_info) ``` ```{code-cell} ipython3 -allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v) -``` - -```{code-cell} ipython3 -present_dict(allocation) -``` - -### round 3 - -```{code-cell} ipython3 -p,bid_info = submit_bid(loser_list, p, ϵ, v, bid_info) -``` - -```{code-cell} ipython3 -p -``` - -```{code-cell} ipython3 -present_dict(bid_info) -``` - -```{code-cell} ipython3 -allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v) +allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v, ϵ) ``` ```{code-cell} ipython3 present_dict(allocation) ``` -### round 4 - -```{code-cell} ipython3 -p,bid_info = submit_bid(loser_list, p, ϵ, v, bid_info) -``` - -```{code-cell} ipython3 -p -``` - -```{code-cell} ipython3 -present_dict(bid_info) -``` - -Notice that Buyer 3 now switches to bid for house 1 having recongized that house 2 is no longer his best option. - -```{code-cell} ipython3 -allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v) -``` - -```{code-cell} ipython3 -present_dict(allocation) -``` +### later rounds -### round 5 +Every remaining round repeats exactly the same two steps, so rather than writing them out one at a time we iterate until no price changes. -```{code-cell} ipython3 -p,bid_info = submit_bid(loser_list, p, ϵ, v, bid_info) -``` +Watch for two things as the rounds unfold: a buyer who is outbid switches to whichever house now offers him the largest surplus, and the auction stops as soon as a round passes in which no losing buyer is willing to raise any price. ```{code-cell} ipython3 -p -``` +round_number = 2 -```{code-cell} ipython3 -present_dict(bid_info) -``` +while True: + round_number += 1 + p_previous = p.copy() -Now Buyer 1 bids for house 1 again with price at 4, which crowds out Buyer 3, marking the end of the auction. + p, bid_info = submit_bid(loser_list, p, ϵ, v, bid_info) + allocation, winner_list, loser_list = check_terminal_condition(bid_info, p, v, ϵ) -```{code-cell} ipython3 -allocation,winner_list,loser_list = check_terminal_condition(bid_info, p, v) -``` + print(f'round {round_number}: prices {p}') + present_dict(allocation) -```{code-cell} ipython3 -present_dict(allocation) + if np.array_equal(p, p_previous): + print('no price changed during this round, so the auction ends') + break ``` ```{code-cell} ipython3 @@ -690,11 +671,13 @@ class ascending_bid_auction: self.v = v.copy() self.n,self.m = self.v.shape - self.r = r + # prices are stored as floats so that a bid increment ϵ < 1 is allowed + self.r = np.asarray(r, dtype=float) self.ϵ = ϵ - self.p = r.copy() + self.p = self.r.copy() self.buyer_list = np.arange(self.m) self.house_list = np.arange(self.n) + self.log = print # replaced by a no-op when verbose=False self.bid_info_history = [] self.allocation_history = [] self.winner_history = [] @@ -726,7 +709,7 @@ class ascending_bid_auction: p_start_mat = (self.ϵ + self.p)[:,None] @ np.ones(self.m)[None,:] self.surplus_value = self.v - p_start_mat - # we only care about active buyers who have positve surplus values + # we only care about active buyers who have positive surplus values active_buyer_diagnosis = (self.surplus_value > 0).any(axis = 0) active_buyer_list = self.buyer_list[active_buyer_diagnosis] active_buyer_surplus_value = self.surplus_value[:,active_buyer_diagnosis] @@ -741,17 +724,17 @@ class ascending_bid_auction: bid_info[house_num] = active_buyer_list[active_buyer_choice == house_num] self.bid_info_history.append(bid_info) - print('The bid information is') + self.log('The bid information is') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *bid_info.keys()] ymtb.add_row(['Buyer', *bid_info.values()]) - print(ymtb) + self.log(ymtb) - print('The bid prices for houses are') + self.log('The bid prices for houses are') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *self.house_list] ymtb.add_row(['Price', *self.p]) - print(ymtb) + self.log(ymtb) self.winner_list=[np.random.choice(bid_info[ii]) for ii in bid_info.keys()] self.winner_history.append(self.winner_list) @@ -763,20 +746,20 @@ class ascending_bid_auction: self.loser_list = list(loser_set) self.loser_history.append(self.loser_list) - print('The winners are') - print(self.winner_list) + self.log('The winners are') + self.log(self.winner_list) - print('The losers are') - print(self.loser_list) - print('\n') + self.log('The losers are') + self.log(self.loser_list) + self.log('\n') def check_terminal_condition(self): loser_num = len(self.loser_list) if loser_num == 0: - print('The auction ends because every buyer gets one house.') - print('\n') + self.log('The auction ends because every buyer gets one house.') + self.log('\n') return True p_mat = (self.ϵ + self.p)[:,None] @ np.ones(loser_num)[None,:] @@ -812,17 +795,17 @@ class ascending_bid_auction: bid_info[house_num] = bid_info_active_loser[house_num] self.bid_info_history.append(bid_info) - print('The bid information is') + self.log('The bid information is') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *bid_info.keys()] ymtb.add_row(['Buyer', *bid_info.values()]) - print(ymtb) + self.log(ymtb) - print('The bid prices for houses are') + self.log('The bid prices for houses are') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *self.house_list] ymtb.add_row(['Price', *self.p]) - print(ymtb) + self.log(ymtb) self.winner_list=[np.random.choice(bid_info[ii]) for ii in bid_info.keys()] self.winner_history.append(self.winner_list) @@ -834,79 +817,81 @@ class ascending_bid_auction: self.loser_list = list(loser_set) self.loser_history.append(self.loser_list) - print('The winners are') - print(self.winner_list) + self.log('The winners are') + self.log(self.winner_list) - print('The losers are') - print(self.loser_list) - print('\n') + self.log('The losers are') + self.log(self.loser_list) + self.log('\n') - def start_auction(self): - print('The Ascending Bid Auction for Houses') - print('\n') + def start_auction(self, verbose=True): + # set verbose=False to run the auction without printing a report + self.log = print if verbose else (lambda *args, **kwargs: None) + self.log('The Ascending Bid Auction for Houses') + self.log('\n') - print('Basic Information: %d houses, %d buyers'%(self.n, self.m)) + self.log('Basic Information: %d houses, %d buyers'%(self.n, self.m)) - print('The valuation matrix is as follows') + self.log('The valuation matrix is as follows') ymtb = pt.PrettyTable() ymtb.field_names = ['Buyer Number', *(np.arange(self.m))] for ii in range(self.n): ymtb.add_row(['House %d'%(ii), *self.v[ii,:]]) - print(ymtb) + self.log(ymtb) - print('The reservation prices for houses are') + self.log('The reservation prices for houses are') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *self.house_list] ymtb.add_row(['Price', *self.r]) - print(ymtb) - print('The minimum increment of bid price is %.2f' % self.ϵ) - print('\n') + self.log(ymtb) + self.log('The minimum increment of bid price is %.2f' % self.ϵ) + self.log('\n') ctr = 1 if self.check_kick_off_condition(): - print('Auction starts successfully') - print('\n') - print('Round %d'% ctr) + self.log('Auction starts successfully') + self.log('\n') + self.log('Round %d'% ctr) self.submit_initial_bid() while True: if self.check_terminal_condition(): - print('Auction ends') - print('\n') + self.log('Auction ends') + self.log('\n') - print('The final result is as follows') - print('\n') - print('The allocation plan is') + self.log('The final result is as follows') + self.log('\n') + self.log('The allocation plan is') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *self.allocation.keys()] ymtb.add_row(['Buyer', *self.allocation.values()]) - print(ymtb) + self.log(ymtb) - print('The bid prices for houses are') + self.log('The bid prices for houses are') ymtb = pt.PrettyTable() ymtb.field_names = ['House Number', *self.house_list] ymtb.add_row(['Price', *self.p]) - print(ymtb) + self.log(ymtb) - print('The winners are') - print(self.winner_list) + self.log('The winners are') + self.log(self.winner_list) - print('The losers are') - print(self.loser_list) + self.log('The losers are') + self.log(self.loser_list) self.house_unsold_list = list(set(self.house_list).difference(set(self.allocation.keys()))) - print('The houses unsold are') - print(self.house_unsold_list) + self.log('The houses unsold are') + self.log(self.house_unsold_list) self.total_revenue = self.p[list(self.allocation.keys())].sum() - print('The total revenue is %.2f' % self.total_revenue) + self.log('The total revenue is %.2f' % self.total_revenue) break ctr += 1 - print('Round %d'% ctr) + self.log('Round %d'% ctr) self.submit_bid() # we compute the surplus matrix S and the quantity matrix X as required in 1.1 @@ -930,7 +915,7 @@ class ascending_bid_auction: self.allocation = dict_temp else: - print('The auction can not start because of high reservation prices') + self.log('The auction can not start because of high reservation prices') ``` Let's use our class to conduct the auction described in one of the above examples. @@ -971,7 +956,7 @@ auction_2 = ascending_bid_auction(v2, r, ϵ) auction_2.start_auction() ``` -**2. multilple excess buyers** +**2. multiple excess buyers** ```{code-cell} ipython3 v3 = np.array([[8,5,9,4,3],[4,11,7,4,6],[9,7,6,4,2]]) @@ -1021,7 +1006,7 @@ auction_6.start_auction() +++ -We now decribe an alternative way for society to allocate $n$ houses to $m$ possible buyers in a way that maximizes +We now describe an alternative way for society to allocate $n$ houses to $m$ possible buyers in a way that maximizes total value across all potential buyers. We continue to assume that each buyer can purchase at most one house. @@ -1030,7 +1015,7 @@ The mechanism is a very special case of a Groves-Clarke mechanism {cite}`Groves Its special structure substantially simplifies writing Python code to find an optimal allocation. -Our mechanims works like this. +Our mechanism works like this. * The values $V_{ij}$ are private information to person $j$ @@ -1040,9 +1025,9 @@ Our mechanims works like this. * The social planner tells no one these, but uses them to allocate houses and set prices -* The mechanism is designed in a way that makes all prospective buyers want to tell the planner their private values +* The mechanism is *intended* to make all prospective buyers want to tell the planner their private values - - truth telling is a dominant strategy for each potential buyer + - we shall put this intention to a test below, and find that it is not achieved * The planner finds a house, bidder pair with highest private value by computing $(\tilde i, \tilde j) = \operatorname{argmax} (V_{ij})$ @@ -1079,7 +1064,7 @@ $$ ```{code-cell} ipython3 np.random.seed(666) -V_orig = np.array([[10, 9, 8, 7, 6], # record the origianl values +V_orig = np.array([[10, 9, 8, 7, 6], # record the original values [9, 9, 7, 6, 6], [8, 6, 6, 9, 4], [7, 5, 6, 4, 9]]) @@ -1235,12 +1220,12 @@ class GC_Mechanism: Find the house-buyer pair with the highest value. When the highest private value corresponds to more than one house, bidder pairs, we choose the pair with the highest sale price. - Moreoever, if the highest sale price corresponds to two or more pairs with highest private value, + Moreover, if the highest sale price corresponds to two or more pairs with highest private value, We randomly choose one. Parameters: ---------- - V: 2 dimensional private value matrix with -1 indicating revomed rows and columns + V: 2 dimensional private value matrix with -1 indicating removed rows and columns Returns: ---------- @@ -1265,7 +1250,12 @@ class GC_Mechanism: return i, j def update_status(self, i, j): - self.p[i] = np.max(np.delete(self.V[i, :], j)) + # the price is the best forgone value of a rival buyer for the same house; + # entries equal to -1 mark buyers and houses already dealt with, and if no + # rival remains there is no competing bid, so the price is zero + rivals = np.delete(self.V[i, :], j) + rivals = rivals[rivals >= 0] + self.p[i] = rivals.max() if rivals.size > 0 else 0.0 self.Q[i, j] = 1 self.V[i, :] = -1 self.V[:, j] = -1 @@ -1273,18 +1263,19 @@ class GC_Mechanism: def calculate_surplus(self): self.S = self.V_orig*self.Q - np.diag(self.p)@self.Q - def start(self): + def start(self, verbose=True): + log = print if verbose else (lambda *args, **kwargs: None) while (np.max(self.V)>=0): i, j = self.find_argmax() self.update_status(i, j) - print("House %i is sold to buyer %i at price %i"%(i[0], j[0], self.p[i[0]])) - print("\n") + log("House %i is sold to buyer %i at price %i"%(i[0], j[0], self.p[i[0]])) + log("\n") self.calculate_surplus() - print("Prices of house:\n", self.p) - print("\n") - print("The status matrix:\n", self.Q) - print("\n") - print("The surplus matrix:\n", self.S) + log("Prices of house:\n", self.p) + log("\n") + log("The status matrix:\n", self.Q) + log("\n") + log("The surplus matrix:\n", self.S) ``` @@ -1337,12 +1328,15 @@ In our setting, equation {eq}`eq:GC2` says that the mechanism charges people for Thus, notice that according to equation {eq}`eq:GC2`: -- unsuccessful prospective buyers pay $0$ because removing them from "society" would not affect the allocation chosen by the mechanim +- unsuccessful prospective buyers pay $0$ because removing them from "society" would not affect the allocation chosen by the mechanism - successful prospective buyers pay the difference between the total value society could achieve without them present and the total value that others present in society do achieve under the mechanism. -The generalized second-price auction described in our pseudo code above does indeed satisfy (1). -We want to compute $\check t_j$ for $j = 1, \ldots, m$ and compare with $p_j$ from the second price auction. +It is natural to ask whether the greedy mechanism that we coded above satisfies {eq}`eq:GC1` and {eq}`eq:GC2`. + +We shall compute $\check t_j$ for $j = 1, \ldots, m$ and compare it with the price $p_j$ that the greedy mechanism actually charges. + +The answer, which we develop in the next two sections, is that it satisfies neither requirement in general. +++ @@ -1410,3 +1404,294 @@ gc_mechanism_exc_2.start() print("\nThe social cost of buyer 2:", np.sum(gc_mechanism_exc_2.Q*gc_mechanism_exc_2.V_orig)-np.sum(np.delete(gc_mechanism.Q*gc_mechanism.V_orig, 2, axis=1))) ``` + +## Is the Greedy Mechanism a VCG Mechanism? + +The mechanism that our benevolent planner runs is **greedy**: it repeatedly grabs the largest remaining entry of the value matrix. + +A genuine VCG mechanism must do two things, namely, satisfy {eq}`eq:GC1` and {eq}`eq:GC2`. + +Let's check them one at a time. + +We'll need a way to compute the value-maximizing allocation, which is a classic **assignment problem** that can be solved by linear programming. + +```{code-cell} ipython3 +from scipy.optimize import linear_sum_assignment + +def optimal_assignment(V): + """ + Solve the assignment problem: allocate at most one house to each buyer + so as to maximize the total value. + Returns the total value and a dictionary mapping buyer -> house. + """ + rows, cols = linear_sum_assignment(-V) + return V[rows, cols].sum(), {int(c): int(r) for r, c in zip(rows, cols)} + +def greedy_outcome(V): + "Run the planner's greedy mechanism quietly and report what it achieves." + gc = GC_Mechanism(V.copy()) + gc.start(verbose=False) + return (gc.Q * gc.V_orig).sum(), gc +``` + +### Does the greedy allocation maximize total value? + +Consider this simple matrix of private values for two houses and two buyers. + +```{code-cell} ipython3 +V_trap = np.array([[10, 9], + [9, 1]]) + +greedy_value, _ = greedy_outcome(V_trap) +best_value, best_alloc = optimal_assignment(V_trap) + +print(f"greedy total value: {greedy_value:.0f}") +print(f"optimal total value: {best_value:.0f} (buyer -> house: {best_alloc})") +``` + +The greedy rule starts from the largest entry, $V_{00} = 10$, and gives house 0 to buyer 0. + +That leaves only $V_{11} = 1$, for a total of $11$. + +But assigning house 0 to buyer 1 and house 1 to buyer 0 yields $9 + 9 = 18$. + +The largest entry was a trap: grabbing it destroyed a better pair of assignments. + +So the greedy mechanism does **not** satisfy {eq}`eq:GC1`. + +{ref}`ha_ex1` asks you to investigate how often this happens. + +The value matrices used earlier in this lecture happen to be ones for which greedy does find the best allocation, which is why the difficulty did not show up before. + +### Are the prices VCG payments? + +Now let's compare the prices that the greedy mechanism charges with the social costs {eq}`eq:GC2`. + +```{code-cell} ipython3 +def vcg_payments(V): + """ + The social cost of each buyer, as defined in {eq}`eq:GC2`: + the total value others could achieve without this buyer, + minus the total value others do achieve when this buyer is present. + """ + total, alloc = optimal_assignment(V) + payments = np.zeros(V.shape[1]) + + for j in range(V.shape[1]): + total_without_j, _ = optimal_assignment(np.delete(V, j, axis=1)) + value_of_others = total - (V[alloc[j], j] if j in alloc else 0.0) + payments[j] = total_without_j - value_of_others + + return payments +``` + +```{code-cell} ipython3 +t_vcg = vcg_payments(V_orig) + +print(f"{'buyer':>6}{'house':>7}{'greedy price':>15}{'VCG payment':>14}") +for j in range(V_orig.shape[1]): + house = [i for i in range(V_orig.shape[0]) if gc_mechanism.Q[i, j] == 1] + if house: + print(f"{j:>6}{house[0]:>7}{gc_mechanism.p[house[0]]:>15.1f}{t_vcg[j]:>14.1f}") +``` + +The prices differ from the VCG payments. + +Notice that the VCG payments are exactly the social costs that we computed in the previous section by excluding one buyer at a time. + +So the planner's mechanism collects **more** than the externality that a winning buyer imposes on others. + +A mechanism that assigns objects greedily and charges each winner the best forgone offer for the same object is known as a **generalized second-price** mechanism. + +Generalized second-price mechanisms are widely used -- they are how search engines sell advertising slots -- and it is well known that they are *not* VCG mechanisms {cite}`EdelmanOstrovskySchwarz2007`. + +### Is truthful reporting a dominant strategy? + +Because the mechanism is not VCG, the argument that makes truth telling a dominant strategy does not apply. + +In fact truth telling is *not* a dominant strategy here, and {ref}`ha_ex2` asks you to demonstrate that by finding a profitable lie. + +### What about the ascending bid auction? + +Here is the striking fact promised earlier. + +Our setting is an assignment market in which each buyer wants at most one house. + +For such markets, {cite}`Leonard1983` and {cite}`DemangeGaleSotomayor1986` showed that the **minimum competitive equilibrium prices** coincide with the payments prescribed by the VCG mechanism. + +Since an ascending auction with straightforward bidding converges to those minimum competitive equilibrium prices, it follows that our ascending bid auction implements the VCG outcome as the bid increment $\epsilon$ becomes small. + +It is the auctioneer, not the benevolent planner, who runs a VCG mechanism here. + +{ref}`ha_ex3` asks you to verify this numerically. + +## Exercises + +```{exercise} +:label: ha_ex1 + +We saw that the greedy mechanism can fail to maximize total value. + +How often does it fail? + +Draw many random value matrices with integer entries, run both the greedy mechanism and `optimal_assignment`, and report + +1. the fraction of matrices for which greedy is suboptimal +1. the largest percentage shortfall that you find + +Then explain in words why a greedy rule goes wrong. +``` + +```{solution-start} ha_ex1 +:class: dropdown +``` + +```{code-cell} ipython3 +rng_ha = np.random.default_rng(0) +trials = 200 +failures, worst = 0, 0.0 + +for _ in range(trials): + n_h, m_b = rng_ha.integers(2, 5), rng_ha.integers(2, 6) + if m_b < n_h: + n_h, m_b = m_b, n_h + V_rand = rng_ha.integers(0, 12, size=(n_h, m_b)) + + greedy_value, _ = greedy_outcome(V_rand) + best_value, _ = optimal_assignment(V_rand) + + if greedy_value < best_value - 1e-9: + failures += 1 + worst = max(worst, (best_value - greedy_value)/best_value) + +print(f"greedy was suboptimal in {failures} of {trials} matrices") +print(f"largest shortfall: {100*worst:.1f} per cent of the attainable total value") +``` + +Greedy fails often, and sometimes badly. + +The reason is that the largest entry of the value matrix need not belong to the best assignment. + +Taking it removes a house *and* a buyer from further consideration, and that can destroy two good matches in order to secure one very good one. + +Choosing the best assignment requires weighing all pairings simultaneously, which is what the linear programming routine does. + +```{solution-end} +``` + +```{exercise} +:label: ha_ex2 + +Show that truthful reporting is not a dominant strategy under the planner's greedy mechanism. + +Use the value matrix `V_orig` of the previous section, and consider buyer 0, whose true values for the three houses are $(10, 9, 8)$. + +1. What surplus does buyer 0 receive when everyone reports truthfully? +1. Search over reports $(a, b, c)$ with integer entries between 0 and 10 that buyer 0 might submit instead, holding fixed the reports of all other buyers, and find the report that maximizes buyer 0's **true** surplus. +1. Explain what the manipulation accomplishes. +``` + +```{solution-start} ha_ex2 +:class: dropdown +``` + +```{code-cell} ipython3 +def surplus_of_buyer_0(V_true, V_reported): + "True surplus of buyer 0 when the mechanism is run on the reported values." + gc = GC_Mechanism(V_reported.copy()) + gc.start(verbose=False) + for i in range(V_true.shape[0]): + if gc.Q[i, 0] == 1: + return V_true[i, 0] - gc.p[i], i, gc.p[i] + return 0.0, None, None # buyer 0 wins nothing + +np.random.seed(0) +truth_surplus, house, price = surplus_of_buyer_0(V_orig, V_orig) +print(f"truthful report: buyer 0 wins house {house} at price {price:.0f}, " + f"surplus {truth_surplus:.0f}") + +best = (truth_surplus, None) +for a in range(11): + for b in range(11): + for c in range(11): + V_lie = V_orig.copy() + V_lie[:, 0] = [a, b, c] + surplus, h, pr = surplus_of_buyer_0(V_orig, V_lie) + if surplus > best[0] + 1e-9: + best = (surplus, (a, b, c, h, pr)) + +report = best[1] +print(f"best misreport found: {report[:3]} wins house {report[3]} " + f"at price {report[4]:.0f}, surplus {best[0]:.0f}") +``` + +By reporting that he places no value on the two expensive houses, buyer 0 steps out of the competition for them. + +The mechanism then awards him a cheaper house, and charges him the best forgone offer for *that* house, which is much less than what he paid when he told the truth. + +His true value for the house he ends up with is high enough that this is a profitable deviation. + +Under a genuine VCG mechanism this could not happen: a buyer's payment never depends on his own report, only on the harm he does to others. + +```{solution-end} +``` + +```{exercise} +:label: ha_ex3 + +Verify that the ascending bid auction implements the VCG outcome. + +Use the value matrix + +$$ +v =\begin{bmatrix} 8 & 5 & 9 & 4 \cr + 4 & 11 & 7 & 4 \cr + 9 & 7 & 6 & 4 \end{bmatrix} +$$ + +with reservation prices of zero, so that prices are free to settle at their competitive levels. + +1. Compute the VCG payments and the value-maximizing allocation. +1. Run the ascending bid auction for $\epsilon = 1, 0.5, 0.25, 0.1, 0.05$ and record what each winning buyer pays. +1. Does the auction maximize total value? Do its payments approach the VCG payments? +``` + +```{solution-start} ha_ex3 +:class: dropdown +``` + +```{code-cell} ipython3 +v_ex = np.array([[8, 5, 9, 4], + [4, 11, 7, 4], + [9, 7, 6, 4]]) + +t_vcg_ex = vcg_payments(v_ex) +best_value_ex, best_alloc_ex = optimal_assignment(v_ex) + +print(f"VCG payments by buyer: {t_vcg_ex}") +print(f"value-maximizing allocation (buyer -> house): {best_alloc_ex}") +print(f"maximum total value: {best_value_ex:.0f}\n") + +print(f"{'ε':>7}{'total value':>14} payments by buyer") +for ϵ_try in (1.0, 0.5, 0.25, 0.1, 0.05): + np.random.seed(0) + auction = ascending_bid_auction(v_ex, np.zeros(3), ϵ_try) + auction.start_auction(verbose=False) + + payments = {j: round(float(auction.p[i]), 3) + for i in range(v_ex.shape[0]) + for j in range(v_ex.shape[1]) if auction.Q[i, j] == 1} + total_value = (auction.Q[:, :v_ex.shape[1]] * v_ex).sum() + print(f"{ϵ_try:>7}{total_value:>14.0f} {payments}") +``` + +The auction always attains the maximum total value. + +As $\epsilon$ shrinks, each winner's payment approaches that buyer's VCG payment. + +This illustrates the theorem of {cite}`Leonard1983` and {cite}`DemangeGaleSotomayor1986`: in an assignment market, the minimum competitive equilibrium prices are exactly the VCG payments. + +The auctioneer never asks anyone to reveal a private value, yet the ascending auction arrives at the outcome that a correctly designed planner would compute from full information. + +```{solution-end} +``` diff --git a/lectures/lq_mean_field_games.md b/lectures/lq_mean_field_games.md new file mode 100644 index 000000000..d017f3e6a --- /dev/null +++ b/lectures/lq_mean_field_games.md @@ -0,0 +1,1644 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: Python 3 + language: python + name: python3 +--- + +(lq_mean_field_games)= +```{raw} jupyter + +``` + +# Linear Quadratic Mean Field Games + +```{contents} Contents +:depth: 2 +``` + +## Overview + +A **mean field game** describes a continuum of small agents, each of whom solves a dynamic optimization problem whose payoff depends on what everybody else is doing, summarized by the cross-sectional distribution of states. + +The framework was introduced by {cite}`LasryLions2007` and, independently, by {cite}`HuangMalhameCaines2006`. + +It pairs two partial differential equations: + +* a **Hamilton-Jacobi-Bellman** equation that runs backward in time and describes an individual's optimal choice, given paths for the aggregates +* a **Kolmogorov forward** equation that runs forward in time and describes how the distribution of individual states evolves, given those choices + +An equilibrium requires that the aggregates that agents take as given are the ones that their own decisions generate. + +The same pair of equations is the workhorse of continuous-time heterogeneous-agent macroeconomics; see {cite}`AchdouEtAl2022`. + +Readers of {doc}`rational_expectations` will recognize that requirement. + +It is the "Big $Y$, little $y$" idea, now applied to an entire distribution rather than to a single number. + +This lecture studies a tractable special case in which the payoff is quadratic and the state evolves linearly, following {cite}`AlvarezArgente2026`. + +Two kinds of interaction appear: + +* agents care about the cross-sectional average **state** $X$, through a matrix $\Theta_X$ +* agents care about the cross-sectional average **action** $\mathcal A$, through a matrix $\Theta_{\mathcal A}$ + +The main result is a striking simplification: + +```{note} +The equilibrium of the mean field game solves the algebraic Riccati equation of a **single-agent** linear quadratic regulator problem, in which the curvature matrices $Q$ and $\Gamma$ are replaced by + +$$ +Q + \Theta_X \qquad\text{and}\qquad \Gamma + \Theta_{\mathcal A} . +$$ +``` + +Everything we know about the linear regulator can therefore be brought to bear on the equilibrium: existence conditions, uniqueness, comparative statics, and numerical methods. + +We then put the framework to work on two economic examples from {cite}`AlvarezArgente2026`: an industry equilibrium with capital accumulation, and a multiproduct pricing problem with Kimball demand. + +Riccati equations appear in several other QuantEcon lectures: + +* {doc}`lqcontrol` introduces the linear regulator and its Riccati equation +* {doc}`lagrangian_lqdp` studies the state-costate system and its stable invariant subspace, which is exactly the structure we meet below +* {doc}`markov_perf` studies dynamic games with *finitely many* players, where each player has a Riccati equation and the equations are coupled +* {doc}`kalman` presents the Riccati equation that is dual to the control problem + +Let's start with some imports: + +```{code-cell} ipython +import numpy as np +import matplotlib.pyplot as plt +from scipy.linalg import solve_continuous_are +``` + +## The environment + +There is a continuum of agents. + +An individual has state $x \in \mathbb R^n$ and takes an action $\alpha \in \mathbb R^k$. + +The cross-sectional average state and action are + +$$ +X = \int x \, m(x) dx , \qquad +\mathcal A = \int \alpha_*(x) \, m(x) dx , +$$ + +where $m$ is the density of the state across agents. + +The period return is quadratic in all four objects: + +```{math} +:label: mfg_return +F(x,X) + R(\alpha, \mathcal A) += -\tfrac12 x' Q x - X'\Theta_X x + -\tfrac12 \alpha' \Gamma \alpha - \mathcal A' \Theta_{\mathcal A}\alpha . +``` + +The individual state follows + +```{math} +:label: mfg_state +dx = (B\alpha - Ax) dt + \Sigma^{1/2} dW , +``` + +where $W$ is an $n$-dimensional Brownian motion and shocks are purely idiosyncratic for now. + +Agents discount at rate $\rho > 0$. + +We assume that $Q$ and $\Gamma$ are positive definite, that $\Sigma$ is positive semi-definite, that $B$ has full row rank, and that $\Theta_X$ and $\Theta_{\mathcal A}$ are symmetric. + +### Complements and substitutes + +From {eq}`mfg_return`, the cross derivative between an agent's own state and the average state is $-\Theta_X$, and between her own action and the average action it is $-\Theta_{\mathcal A}$. + +So + +* states are **strategic complements** when $-\Theta_X$ is positive definite, i.e. when $\Theta_X$ is negative definite +* states are **strategic substitutes** when $\Theta_X$ is positive definite + +and similarly for actions. + +In the Loewner ordering, *smaller* $\Theta$ means *more* complementarity. + +```{note} +The monotonicity condition that {cite}`LasryLions2007` impose to obtain uniqueness corresponds here to $-\Theta_X$ being negative semi-definite, that is, to strategic *substitutability* in states. + +We will not need it: in this linear quadratic setting there is at most one equilibrium whether interactions are complements or substitutes. +``` + +### Equilibrium + +Taking the paths $\{X(t), \mathcal A(t)\}$ as given, an agent's value function satisfies the HJB equation + +```{math} +:label: mfg_hjb +\rho u(x,t) = -\tfrac12 x'Qx - X(t)'\Theta_X x + + H(u_x(x,t), x, \mathcal A(t)) + + \tfrac12 \operatorname{tr}(\Sigma u_{xx}(x,t)) + u_t(x,t) , +``` + +where the Hamiltonian is + +$$ +H(p, x, \mathcal A) = \max_{\alpha} +\left\{ -\tfrac12\alpha'\Gamma\alpha - \mathcal A'\Theta_{\mathcal A}\alpha ++ p'(B\alpha - Ax) \right\} , +$$ + +with maximizer $\alpha_*(p,\mathcal A) = \Gamma^{-1}(B'p - \Theta_{\mathcal A}\mathcal A)$. + +The density evolves according to the Kolmogorov forward equation + +```{math} +:label: mfg_kfe +m_t(x,t) = -\operatorname{div}\left( H_p(u_x(x,t),x,\mathcal A(t)) m(x,t)\right) ++ \tfrac12 \operatorname{tr}(\Sigma m_{xx}(x,t)) , +``` + +and an **equilibrium** is a value function, a density, and paths for $X$ and $\mathcal A$ that satisfy {eq}`mfg_hjb`, {eq}`mfg_kfe`, and the consistency requirements that $X$ and $\mathcal A$ really are the cross-sectional averages implied by $m$ and by the optimal policy. + +## Solving the individual problem + +Given the aggregate paths, an individual faces a time-varying linear quadratic regulator problem, so her value function is quadratic: + +$$ +u(x,t) = \beta_0(t) + \beta_1(t)'x + \tfrac12 x'\beta_2(t)x . +$$ + +Substituting into {eq}`mfg_hjb` and matching terms of each order gives three differential equations. + +The one for $\beta_2$ is + +```{math} +:label: mfg_beta2_ode +\dot\beta_2 = Q - \beta_2 B\Gamma^{-1}B'\beta_2 + \beta_2 A + A'\beta_2 + \rho\beta_2 . +``` + +Notice what is *absent* from {eq}`mfg_beta2_ode`: neither interaction matrix appears. + +The curvature of an individual's value function is therefore the same as it would be if she were alone in the world. + +Because the individual problem is concave and stationary, $\beta_2(t)$ equals the constant $\bar\beta_2$, the negative definite solution of + +```{math} +:label: mfg_beta2 +\bar\beta_2 B\Gamma^{-1}B'\bar\beta_2 = Q + \rho\bar\beta_2 + \bar\beta_2 A + A'\bar\beta_2 . +``` + +This is the familiar algebraic Riccati equation of {doc}`lqcontrol`, written in continuous time. + +The optimal action is + +$$ +\alpha_*(x,t) = \Gamma^{-1}\left[B'(\beta_1(t) + \bar\beta_2 x) - \Theta_{\mathcal A}\mathcal A(t)\right] . +$$ + +Averaging across agents and solving the resulting fixed point in $\mathcal A$ gives + +```{math} +:label: mfg_aggregate_action +\mathcal A(t) = (\Gamma + \Theta_{\mathcal A})^{-1} B' \left(\beta_1(t) + \bar\beta_2 X(t)\right) . +``` + +Equation {eq}`mfg_aggregate_action` is where the action interaction first bites: each agent responds to the average action, and solving for the average that is consistent with everyone doing so replaces $\Gamma$ by $\Gamma + \Theta_{\mathcal A}$. + +## A state-costate system + +Two objects now remain: the linear coefficient $\beta_1(t)$ of the value function, and the aggregate state $X(t)$. + +Differentiating and aggregating gives a pair of linear differential equations, + +```{math} +:label: mfg_hamiltonian_system +\begin{bmatrix} \dot\beta_1 \\ \dot X \end{bmatrix} += \mathcal H \begin{bmatrix} \beta_1 \\ X\end{bmatrix}, +\qquad +\mathcal H = +\begin{bmatrix} +\rho I + A' - \bar\beta_2 \Lambda & \Theta_X + \bar\beta_2(B\Gamma^{-1}B' - \Lambda)\bar\beta_2 \\ +\Lambda & -A + \Lambda\bar\beta_2 +\end{bmatrix}, +``` + +where we abbreviate + +$$ +\Lambda \equiv B(\Gamma + \Theta_{\mathcal A})^{-1}B' . +$$ + +This is a **state-costate** system of exactly the kind studied in {doc}`lagrangian_lqdp`. + +The aggregate state $X$ has an initial condition, namely the mean of the initial distribution. + +The costate $\beta_1$ does not: it must be chosen so that the solution does not violate the agent's transversality condition, which here requires the eigenvalues governing the path to have real parts below $\rho/2$. + +Because $\Theta_X$, $B\Gamma^{-1}B'$ and $\Lambda$ are symmetric, $\mathcal H - \tfrac\rho2 I$ is a Hamiltonian matrix, so its eigenvalues are symmetric about the origin. + +Equivalently: + +```{prf:proposition} +:label: mfg_prop_roots + +If $\lambda$ is an eigenvalue of $\mathcal H$, then so are $\rho - \lambda$, $\bar\lambda$, and $\rho - \bar\lambda$. +``` + +Exactly $n$ eigenvalues can therefore have real parts below $\rho/2$, which pins down a unique stable invariant subspace and hence at most one equilibrium. + +This is the standard connection between algebraic Riccati equations and invariant subspaces of Hamiltonian matrices, treated at length by {cite}`LancasterRodman1995`. + +## The equilibrium Riccati equation + +We look for a saddle path along which the costate is a linear function of the state, $\beta_1(t) = S X(t)$. + +Substituting into {eq}`mfg_hamiltonian_system` gives a quadratic matrix equation for $S$ whose coefficients involve $\bar\beta_2$. + +That equation looks forbidding, but a change of variable transforms it. + +Define + +$$ +P \equiv S + \bar\beta_2 . +$$ + +```{prf:proposition} +:label: mfg_prop_riccati + +An equilibrium is characterized by a matrix $P$ solving + +$$ +P \, B(\Gamma+\Theta_{\mathcal A})^{-1}B' \, P += Q + \Theta_X + \rho P + PA + A'P , +$$ + +with aggregate dynamics + +$$ +\dot X = \left(B(\Gamma+\Theta_{\mathcal A})^{-1}B'P - A\right) X . +$$ + +The equilibrium requires all eigenvalues of the closed-loop matrix to have real parts below $\rho/2$. +``` + +Compare this with the single-agent equation {eq}`mfg_beta2`. + +They have **the same form**. + +The only difference is that $Q$ has become $Q + \Theta_X$ and $\Gamma$ has become $\Gamma + \Theta_{\mathcal A}$. + +```{prf:proposition} +:label: mfg_prop_equivalence + +The equilibrium of a linear quadratic mean field game with interaction matrices $\Theta_X$ and $\Theta_{\mathcal A}$, and its aggregate law of motion, coincide with the solution of a single-agent linear quadratic control problem whose state curvature is $Q+\Theta_X$ and whose action curvature is $\Gamma+\Theta_{\mathcal A}$. +``` + +This is the organizing result of the lecture. + +It says that strategic interaction does not change the *form* of the problem that determines aggregate dynamics; it changes the *curvatures* that enter it. + +An immediate economic implication is that two models with very different strategic interactions can generate identical aggregate dynamics, provided the effective curvatures agree. + +### Existence and uniqueness + +Since $P$ solves a standard Riccati equation, standard conditions apply. + +Define + +```{math} +:label: mfg_E +E \equiv Q + \Theta_X + \left(A' + \tfrac\rho2 I\right) +\left[B(\Gamma+\Theta_{\mathcal A})^{-1}B'\right]^{-1} +\left(A + \tfrac\rho2 I\right) . +``` + +```{prf:proposition} +:label: mfg_prop_existence + +Suppose $B(\Gamma+\Theta_{\mathcal A})^{-1}B'$ is invertible. + +1. A necessary condition for an equilibrium is that $E$ be positive semi-definite. +1. If $Q + \Theta_X$ and $\Gamma + \Theta_{\mathcal A}$ are positive definite, an equilibrium exists and is unique. +1. There is at most one equilibrium. +``` + +The sufficient condition has a clean reading: **effective** curvature must remain positive in both states and actions. + +Complementarity is therefore permissible, but only up to a point. + +Notice also that the two interactions enter $E$ additively, so substitutability in actions can offset complementarity in states. + +## Computing equilibria + +`scipy.linalg.solve_continuous_are(A, B, Q, R)` returns the stabilizing solution $X$ of + +$$ +A'X + XA - XBR^{-1}B'X + Q = 0 . +$$ + +Our equation differs in two ways: our $P$ is negative definite, and we discount. + +Writing $P = -\tilde P$ and collecting terms shows that our equation is the standard one with $A$ replaced by $-(A + \tfrac\rho2 I)$. + +The discount rate enters exactly as it does in the "$\rho/2$ shift" familiar from continuous-time control. + +```{code-cell} ipython +def mfg_riccati(A, B, Q, Γ, ρ): + """ + Solve P B Γ^{-1} B' P = Q + ρ P + P A + A' P for the negative + definite stabilizing solution P. + + Passing Q + Θ_X and Γ + Θ_A gives the equilibrium of the mean field game; + passing Q and Γ gives the single agent's value function curvature. + """ + n = A.shape[0] + return -solve_continuous_are(-(A + ρ/2 * np.eye(n)), B, Q, Γ) + +def closed_loop(A, B, Γ_eff, P): + "The matrix governing aggregate dynamics, Ẋ = (B Γ_eff^{-1} B' P - A) X." + return B @ np.linalg.solve(Γ_eff, B.T) @ P - A +``` + +Let's set up a two-dimensional example with complementarity in states and substitutability in actions, the configuration that {cite}`AlvarezArgente2026` obtain from an industry equilibrium with capital accumulation. + +```{code-cell} ipython +ρ = 0.05 +A = np.array([[0.4, 0.1], + [0.0, 0.3]]) +B = np.eye(2) +Q = np.array([[1.0, 0.2], + [0.2, 0.8]]) +Γ = np.array([[1.0, 0.1], + [0.1, 1.2]]) +Θ_X = np.array([[-0.3, 0.05], # negative definite: complements in states + [0.05, -0.2]]) +Θ_A = np.array([[0.2, 0.0], # positive definite: substitutes in actions + [0.0, 0.1]]) + +β2 = mfg_riccati(A, B, Q, Γ, ρ) # single agent +P = mfg_riccati(A, B, Q + Θ_X, Γ + Θ_A, ρ) # equilibrium + +print("individual curvature β̄₂ =\n", β2.round(4)) +print("\nequilibrium matrix P =\n", P.round(4)) +``` + +Let's verify that $P$ really does solve the equilibrium Riccati equation, and compare the aggregate dynamics with what a lone agent would choose. + +```{code-cell} ipython +Λ = B @ np.linalg.solve(Γ + Θ_A, B.T) +residual = P @ Λ @ P - (Q + Θ_X + ρ*P + P @ A + A.T @ P) +print(f"Riccati residual: {np.abs(residual).max():.2e}") + +G = closed_loop(A, B, Γ, β2) # dynamics without any interaction +JG = closed_loop(A, B, Γ + Θ_A, P) # equilibrium dynamics + +print("\neigenvalues without interactions:", np.linalg.eigvals(G).round(4)) +print("eigenvalues in equilibrium: ", np.linalg.eigvals(JG).round(4)) +``` + +The equilibrium eigenvalues are closer to zero, so aggregate adjustment is slower than it would be if each agent ignored everyone else. + +### Checking the Hamiltonian structure + +{prf:ref}`mfg_prop_roots` says the eigenvalues of $\mathcal H$ come in pairs $\{\lambda, \rho - \lambda\}$, and that the $n$ eigenvalues with real parts below $\rho/2$ are the ones that govern equilibrium dynamics. + +Let's check both claims. + +```{code-cell} ipython +BΓB = B @ np.linalg.solve(Γ, B.T) +n = A.shape[0] + +H = np.block([[ρ*np.eye(n) + A.T - β2 @ Λ, Θ_X + β2 @ (BΓB - Λ) @ β2], + [Λ, -A + Λ @ β2]]) + +ev = np.linalg.eigvals(H) +print("eigenvalues of ℋ:", np.sort(ev.real).round(4)) +print("paired as λ and ρ - λ:", + np.allclose(np.sort(ev.real), np.sort(ρ - ev.real))) + +stable = np.sort(ev.real[ev.real < ρ/2]) +print("\nstable half of ℋ: ", stable.round(4)) +print("closed-loop eigenvalues:", np.sort(np.linalg.eigvals(JG).real).round(4)) +``` + +The saddle path is the stable invariant subspace of $\mathcal H$, exactly as in {doc}`lagrangian_lqdp`. + +## The scalar case + +With $n = k = 1$ everything is explicit. + +Write $q, a, b, \gamma, \theta_X, \theta_{\mathcal A}$ for the scalars. + +The Riccati equation becomes a quadratic, and the admissible root gives the aggregate eigenvalue + +```{math} +:label: mfg_scalar_lambda +\lambda = \frac\rho2 - \sqrt{\left(\frac\rho2 + a\right)^2 ++ \frac{b^2(q + \theta_X)}{\gamma + \theta_{\mathcal A}}} . +``` + +An equilibrium exists if and only if the term under the square root is positive, + +```{math} +:label: mfg_scalar_existence +q + \theta_X + \left(\frac\rho2+a\right)^2\frac{\gamma+\theta_{\mathcal A}}{b^2} > 0 , +``` + +and the equilibrium is stable, meaning $\lambda<0$, if and only if + +```{math} +:label: mfg_scalar_stability +q + \theta_X + a(a+\rho)\frac{\gamma+\theta_{\mathcal A}}{b^2} > 0 . +``` + +Formula {eq}`mfg_scalar_lambda` displays the two comparative statics at a glance. + +More complementarity in **states** (a smaller $\theta_X$) raises $\lambda$ and makes aggregate dynamics **more** persistent: when others stay away from the steady state, each agent has less reason to return to it. + +More complementarity in **actions** (a smaller $\theta_{\mathcal A}$) lowers $\lambda$ and makes dynamics **less** persistent: when others adjust, each agent wants to adjust too. + +Let's confirm that our solver reproduces {eq}`mfg_scalar_lambda`. + +```{code-cell} ipython +def scalar_lambda(ρ, a, b, q, γ, θ_X, θ_A): + "Closed-form aggregate eigenvalue in the scalar case." + return ρ/2 - np.sqrt((ρ/2 + a)**2 + b**2*(q + θ_X)/(γ + θ_A)) + +ρ_s, a, b, q, γ = 0.05, 0.3, 1.0, 1.0, 1.0 + +print(f"{'θ_X':>6}{'θ_A':>6}{'solver':>12}{'closed form':>14}") +for θ_X, θ_A in ((0.0, 0.0), (-0.5, 0.0), (0.0, 0.5), (-0.5, 0.5)): + P_s = mfg_riccati(np.array([[a]]), np.array([[b]]), + np.array([[q + θ_X]]), np.array([[γ + θ_A]]), ρ_s) + λ_num = closed_loop(np.array([[a]]), np.array([[b]]), + np.array([[γ + θ_A]]), P_s)[0, 0] + print(f"{θ_X:>6}{θ_A:>6}{λ_num:>12.6f}{scalar_lambda(ρ_s,a,b,q,γ,θ_X,θ_A):>14.6f}") +``` + +## An industry equilibrium with capital accumulation + +The scalar model is not a toy. + +{cite}`AlvarezArgente2026` show that it describes an industry equilibrium with capital accumulation in which **both** kinds of interaction appear, with opposite signs. + +There is a continuum of monopolistically competitive firms. + +A firm with capital $k$ produces $y = k^\nu$ with $0 < \nu < 1$, and a constant returns sector aggregates the differentiated goods with elasticity of substitution $\eta > 1$. + +With the final good as numeraire, a firm that produces $y$ when industry output is $Y$ earns revenue proportional to $y^{1-1/\eta}Y^{1/\eta}$. + +So if all other firms hold capital $K$, operating profit is proportional to + +```{math} +:label: mfg_capital_profit +\Pi(k,K) = k^{\nu(1 - 1/\eta)} K^{\nu/\eta} . +``` + +Capital evolves according to + +$$ +dk = (i - \delta k)dt + k \sigma dW , +$$ + +the firm buys investment goods at price $\mathcal P(I)$, where $I$ is aggregate investment, and it pays a convex adjustment cost $\psi(i)$. + +Let $\bar i = \delta \bar k$ be steady-state investment, normalize $\mathcal P(\bar i) = 1$ and $\psi'(\bar i) = 0$, and write percentage deviations from the deterministic steady state as + +$$ +x = \frac{k - \bar k}{\bar k}, \qquad +X = \frac{K - \bar k}{\bar k}, \qquad +\alpha = \frac{i - \bar i}{\bar i}, \qquad +\mathcal A = \frac{I - \bar i}{\bar i} . +$$ + +A second-order expansion of the return around the steady state, with the objective normalized by $\bar\Pi \equiv \Pi(\bar k, \bar k)$, delivers exactly {eq}`mfg_return` with + +$$ +q = -\frac{\bar k^2 \Pi_{kk}}{\bar\Pi}, \qquad +\theta_X = -\frac{\bar k^2 \Pi_{kK}}{\bar\Pi}, \qquad +\gamma = \frac{\delta^2\bar k^2 \psi''(\bar i)}{\bar\Pi}, \qquad +\theta_{\mathcal A} = \frac{\delta^2\bar k^2 \mathcal P'(\bar i)}{\bar\Pi} , +$$ + +while the state equation becomes $dx = \delta(\alpha - x)dt + \sigma dW$ to first order, so that + +$$ +a = b = \delta . +$$ + +Differentiating {eq}`mfg_capital_profit` gives closed forms: + +```{math} +:label: mfg_capital_coeffs +\begin{aligned} +q &= \nu\frac{\eta-1}{\eta^2}\left[\eta(1-\nu) + \nu\right] > 0 , \\ +\theta_X &= -\frac{\eta-1}{\eta}\frac{\nu^2}{\eta} < 0 , \\ +q + \theta_X &= \frac{\eta-1}{\eta}\nu(1-\nu) > 0 . +\end{aligned} +``` + +Three features of {eq}`mfg_capital_coeffs` deserve emphasis. + +First, $\theta_X < 0$: capital stocks are strategic **complements**, because a larger industry capital stock raises the demand shifter $Y^{1/\eta}$ and hence the marginal profitability of a firm's own capital. + +Second, if the supply of investment goods slopes upward then $\mathcal P'(\bar i) > 0$ and so $\theta_{\mathcal A} > 0$: investment rates are strategic **substitutes**, because everyone investing at once bids up the price of capital goods. + +Third, $q + \theta_X > 0$ for every $\eta > 1$ and every $\nu \in (0,1)$, so by {prf:ref}`mfg_prop_existence` an equilibrium exists and is unique no matter how strong market power is and no matter how close returns to scale come to constant. + +Let's put these formulas into code, and check them against numerical derivatives of the profit function itself. + +```{code-cell} ipython +def cap_q(η, ν): + "Own-state curvature in the capital accumulation example." + return ν*(η - 1)/η**2*(η*(1 - ν) + ν) + +def cap_θ_X(η, ν): + "State interaction in the capital accumulation example." + return -(η - 1)/η*ν**2/η + +η_c, ν_c = 4.0, 0.7 + +Π = lambda k, K: k**(ν_c*(1 - 1/η_c))*K**(ν_c/η_c) + +h = 1e-5 +Π_kk = (Π(1+h, 1) - 2*Π(1, 1) + Π(1-h, 1))/h**2 +Π_kK = (Π(1+h, 1+h) - Π(1+h, 1-h) - Π(1-h, 1+h) + Π(1-h, 1-h))/(4*h**2) + +print(f"{'':>10}{'finite difference':>20}{'closed form':>15}") +print(f"{'q':>10}{-Π_kk/Π(1, 1):>20.6f}{cap_q(η_c, ν_c):>15.6f}") +print(f"{'θ_X':>10}{-Π_kK/Π(1, 1):>20.6f}{cap_θ_X(η_c, ν_c):>15.6f}") +print(f"{'q + θ_X':>10}{-(Π_kk + Π_kK)/Π(1, 1):>20.6f}" + f"{(η_c - 1)/η_c*ν_c*(1 - ν_c):>15.6f}") +``` + +### Calibration + +Take annual units, $\delta = 0.10$, $\rho = 0.05$, an elasticity of substitution $\eta = 4$, and returns to scale $\nu = 0.7$. + +For the two curvatures that the technology does not pin down we set $\gamma = 0.05$, which makes a firm that ignores the industry close half of a capital gap in three years, and $\theta_{\mathcal A} = 0.09$. + +{ref}`mfg_ex5` derives both numbers from an adjustment cost function and an investment supply curve. + +```{code-cell} ipython +ρ_c, δ_c = 0.05, 0.10 +γ_c, θ_A_c = 0.05, 0.09 + +q_c, θ_X_c = cap_q(η_c, ν_c), cap_θ_X(η_c, ν_c) + +def half_life(λ): + "Time for the aggregate state to close half of a gap." + return np.log(2)/(-λ) + +λ_c = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_c, θ_X_c, θ_A_c) + +P_c = mfg_riccati(np.array([[δ_c]]), np.array([[δ_c]]), + np.array([[q_c + θ_X_c]]), np.array([[γ_c + θ_A_c]]), ρ_c) +λ_c_solver = closed_loop(np.array([[δ_c]]), np.array([[δ_c]]), + np.array([[γ_c + θ_A_c]]), P_c)[0, 0] + +print(f"q = {q_c:.4f}, θ_X = {θ_X_c:.4f}, q + θ_X = {q_c + θ_X_c:.4f}") +print(f"λ from the solver = {λ_c_solver:.6f}") +print(f"λ from the closed form = {λ_c:.6f}") +print(f"half-life of aggregate capital = {half_life(λ_c):.2f} years") +``` + +### How much do the interactions matter? + +Because {eq}`mfg_scalar_lambda` depends on $\theta_X$ and $\theta_{\mathcal A}$ separately, we can switch each interaction off and read the answer. + +```{code-cell} ipython +cases = {'no interactions': (0.0, 0.0), + 'state complementarity only': (θ_X_c, 0.0), + 'action substitutability only': (0.0, θ_A_c), + 'equilibrium': (θ_X_c, θ_A_c), + 'planner (both doubled)': (2*θ_X_c, 2*θ_A_c)} + +print(f"{'':>30}{'λ':>10}{'half-life':>12}") +for label, (tx, ta) in cases.items(): + λ_case = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_c, tx, ta) + print(f"{label:>30}{λ_case:>10.4f}{half_life(λ_case):>12.2f}") +``` + +Both interactions slow aggregate capital adjustment, and they do so for different reasons. + +Complementarity in states means that a firm has less reason to rebuild its capital while the rest of the industry is still below its steady state. + +Substitutability in actions means that a burst of industry investment is expensive, so firms spread their investment over time. + +Together they stretch the half-life of the industry's capital stock from three years to five. + +The planner, who internalizes both externalities, is slower still. + +### Micro and macro adjustment speeds + +A distinctive prediction of the model is that individual capital reverts to the mean faster than aggregate capital does, because {eq}`mfg_beta2` contains no interaction matrix. + +```{code-cell} ipython +β2_c = mfg_riccati(np.array([[δ_c]]), np.array([[δ_c]]), + np.array([[q_c]]), np.array([[γ_c]]), ρ_c) +λ_micro = closed_loop(np.array([[δ_c]]), np.array([[δ_c]]), + np.array([[γ_c]]), β2_c)[0, 0] + +print(f"individual half-life = {half_life(λ_micro):.2f} years") +print(f"aggregate half-life = {half_life(λ_c):.2f} years") +print(f"ratio = {half_life(λ_c)/half_life(λ_micro):.2f}") +``` + +A researcher who estimated the speed of capital adjustment from firm-level data, and then used it to predict how quickly the industry responds to an industry-wide shock, would be too fast by two thirds. + +The gap is a pure interaction effect: the same technology and the same adjustment costs generate both numbers. + +### Market power and returns to scale + +How does the industry's adjustment speed depend on the two technological parameters? + +Both enter only through $q + \theta_X = \frac{\eta-1}{\eta}\nu(1-\nu)$, which rises with $\eta$ and is maximized at $\nu = 1/2$. + +```{code-cell} ipython +η_grid = np.linspace(1.2, 12, 200) +ν_grid = np.linspace(0.05, 0.995, 200) + +def macro_micro(η, ν): + "Aggregate and individual half-lives as functions of (η, ν)." + q, θ_X = cap_q(η, ν), cap_θ_X(η, ν) + λ_agg = scalar_lambda(ρ_c, δ_c, δ_c, q, γ_c, θ_X, θ_A_c) + λ_ind = scalar_lambda(ρ_c, δ_c, δ_c, q, γ_c, 0.0, 0.0) + return half_life(λ_agg), half_life(λ_ind) + +fig, axes = plt.subplots(1, 2, figsize=(11, 4.2)) + +hl_agg, hl_ind = np.array([macro_micro(η, ν_c) for η in η_grid]).T +axes[0].plot(η_grid, hl_agg, lw=2, label='industry') +axes[0].plot(η_grid, hl_ind, lw=2, ls='--', label='single firm') +axes[0].set_xlabel('elasticity of substitution $\\eta$') + +hl_agg, hl_ind = np.array([macro_micro(η_c, ν) for ν in ν_grid]).T +axes[1].plot(ν_grid, hl_agg, lw=2, label='industry') +axes[1].plot(ν_grid, hl_ind, lw=2, ls='--', label='single firm') +axes[1].axhline(np.log(2)/δ_c, color='k', lw=1, alpha=0.6) +axes[1].set_xlabel('returns to scale $\\nu$') +axes[1].annotate('$\\ln 2/\\delta$', (0.12, np.log(2)/δ_c - 0.45)) + +for ax in axes: + ax.set_ylabel('half-life in years') + ax.legend() +plt.tight_layout() +plt.show() +``` + +More market power, meaning a smaller $\eta$, makes the industry slower: it weakens the own-capital curvature relative to the interaction, and the left panel shows the industry half-life rising as $\eta$ falls. + +The right panel contains a sharper result. + +As $\nu \to 1$ the effective curvature $q + \theta_X$ vanishes, and {eq}`mfg_scalar_lambda` collapses to $\lambda \to -\delta$: the industry's capital stock then returns to its steady state only through depreciation, with aggregate investment not responding at all. + +{ref}`mfg_ex6` asks you to verify this limit and to explain it. + +### Where are we in the four regions? + +{ref}`mfg_ex1` divides the scalar model into four regions using the thresholds $-\theta^{*}$ and $-\theta^{**}$. + +Since the technology delivers $q + \theta_X > 0$ and an upward sloping investment supply delivers $\theta_{\mathcal A} > 0$, this example is always in the first region. + +The thresholds show how much room there is to spare. + +```{code-cell} ipython +mθ_star = q_c + δ_c*(δ_c + ρ_c)*(γ_c + θ_A_c)/δ_c**2 +mθ_2star = q_c + (ρ_c/2 + δ_c)**2*(γ_c + θ_A_c)/δ_c**2 + +print(f"complementarity in the calibration, -θ_X = {-θ_X_c:.4f}") +print(f"instability threshold, -θ* = {mθ_star:.4f}") +print(f"nonexistence threshold, -θ** = {mθ_2star:.4f}") +``` + +Complementarity would have to be five times stronger than the technology implies before the industry's capital stock stopped converging. + +### Transition paths + +Finally, let's trace out the industry's response to a capital stock that starts ten percent above its steady state. + +Aggregate investment follows from {eq}`mfg_aggregate_action`, which in the scalar case gives $\mathcal A(t) = \frac{\lambda + \delta}{\delta}X(t)$. + +```{code-cell} ipython +t_c = np.linspace(0, 25, 300) +X0_c = 0.10 + +fig, axes = plt.subplots(1, 2, figsize=(11, 4.2)) +for label in ('no interactions', 'equilibrium', 'planner (both doubled)'): + tx, ta = cases[label] + λ_case = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_c, tx, ta) + path = X0_c*np.exp(λ_case*t_c) + axes[0].plot(t_c, 100*path, lw=2, label=label) + axes[1].plot(t_c, 100*(λ_case + δ_c)/δ_c*path, lw=2, label=label) + +axes[0].set_ylabel('capital, % above steady state') +axes[1].set_ylabel('investment, % above steady state') +for ax in axes: + ax.set_xlabel('years') + ax.axhline(0, color='k', lw=0.8) +axes[0].legend() +plt.tight_layout() +plt.show() +``` + +The right panel is the industry's investment response, and it is where the two interactions show up most clearly. + +A firm that ignored the industry would cut investment on impact by thirteen percent of its steady-state level; in equilibrium the cut is under four percent, and it is undone much more slowly. + +The costate $\beta_1(t) + \bar\beta_2 x$ that supports these paths is the marginal value of installed capital, that is, Tobin's $q$. + +{doc}`optimal_growth_uncertainty` studies that object in a setting where the constraint $i \geq 0$ binds and the value function is not differentiable everywhere, which is exactly the nonlinearity that the quadratic approximation here sets aside. + + +## Persistence + +How do the interactions change aggregate adjustment when $n > 1$? + +Stronger complementarity of either kind makes the stabilizing solution $P$ less negative definite. + +That ordering is not enough to sign the change in every eigenvalue of the closed-loop matrix, but it does sign the change in their **sum**. + +Stronger complementarity in states raises the trace of the closed-loop matrix, while stronger complementarity in actions lowers it. + +The trace measures the rate at which a set of initial aggregate states contracts in volume as each point follows its equilibrium path, so in a stable equilibrium the first force slows the collapse toward the steady state and the second speeds it up. + +```{code-cell} ipython +print(f"{'scale on Θ_X':>14}{'trace':>10} eigenvalues") +for scale in (0.0, 0.5, 1.0, 1.5): + P_s = mfg_riccati(A, B, Q + scale*Θ_X, Γ + Θ_A, ρ) + JG_s = closed_loop(A, B, Γ + Θ_A, P_s) + print(f"{scale:>14}{np.trace(JG_s):>10.4f} {np.linalg.eigvals(JG_s).real.round(4)}") +``` + +Raising the scale makes $\Theta_X$ more negative, that is, complementarity stronger, and the trace rises as claimed. + +{ref}`mfg_ex2` asks you to investigate whether *each* eigenvalue must move in the same direction. + +## The planner + +A utilitarian planner internalizes the effect of each agent's state and action on the corresponding cross-sectional averages. + +When $\Theta_X$ and $\Theta_{\mathcal A}$ are symmetric, differentiating the planner's objective doubles each interaction term. + +```{prf:proposition} +:label: mfg_prop_planner + +The planner's allocation coincides with the decentralized equilibrium of an economy whose interaction matrices are $2\Theta_X$ and $2\Theta_{\mathcal A}$, so the planner's Riccati equation is + +$$ +P^{*} B(\Gamma + 2\Theta_{\mathcal A})^{-1}B'P^{*} += Q + 2\Theta_X + \rho P^{*} + P^{*}A + A'P^{*} . +$$ +``` + +Internalizing state complementarity makes the allocation *more* persistent, while internalizing action complementarity makes it *less* persistent. + +When both are present the comparison is ambiguous, as {ref}`mfg_ex3` explores. + +```{code-cell} ipython +P_planner = mfg_riccati(A, B, Q + 2*Θ_X, Γ + 2*Θ_A, ρ) +JG_planner = closed_loop(A, B, Γ + 2*Θ_A, P_planner) + +print("equilibrium eigenvalues:", np.linalg.eigvals(JG).real.round(4)) +print("planner eigenvalues: ", np.linalg.eigvals(JG_planner).real.round(4)) +``` + +Let's see what this means for the path of the aggregate state. + +```{code-cell} ipython +from scipy.linalg import expm + +X0 = np.array([1.0, 0.5]) +times = np.linspace(0, 12, 200) + +paths = {'no interactions': G, 'equilibrium': JG, 'planner': JG_planner} + +fig, ax = plt.subplots(figsize=(8, 4.5)) +for label, M in paths.items(): + traj = np.array([expm(M*t) @ X0 for t in times]) + ax.plot(times, traj[:, 0], lw=2, label=label) +ax.set_xlabel('$t$') +ax.set_ylabel('first component of $X(t)$') +ax.legend() +plt.tight_layout() +plt.show() +``` + +## Multiproduct price setting with Kimball demand + +Our second example is multidimensional, and it delivers a surprise. + +{cite}`AlvarezArgente2026` study a continuum of "stores", each selling $n$ products with constant marginal costs $z_j$. + +Products within a store are aggregated by a CES price index with elasticity $\bar\eta_d$, and stores are aggregated by a symmetric Kimball aggregator {cite}`Kimball1995`. + +Let $\eta_D(y)$ be the elasticity of a store's demand with respect to its relative price, and let $\bar\eta_D > 1$ and $\bar\eta_D'$ denote its level and its derivative at the symmetric point. + +The derivative $\bar\eta_D'$ is the **superelasticity** of demand, and it is the reason Kimball demand is so widely used in models of price setting: a positive superelasticity means that a store which raises its prices above the average faces a more elastic demand, which discourages it from moving away from the crowd. + +That is strategic complementarity in prices, and it is the mechanism that {cite}`KlenowWillis2016` and much of the subsequent literature rely on for real rigidity. + +Writing $x$ and $X$ for log deviations of a store's prices and of the average store's prices from the flexible-price level $\bar p_i = \bar z_i \bar\eta_D/(\bar\eta_D-1)$, and $\bar s$ for the vector of steady-state expenditure shares, the curvature matrices are + +```{math} +:label: mfg_kimball_Q +\begin{aligned} +Q &= (\bar\eta_D - 1)\left[\left(\bar\eta_D - \bar\eta_d ++ \frac{\bar\eta_D'}{\bar\eta_D-1}\right)\bar s\bar s' ++ \bar\eta_d \operatorname{diag}(\bar s)\right] , \\ +\Theta_X &= -\bar\eta_D' \, \bar s \bar s' , \\ +Q + \Theta_X &= (\bar\eta_D-1)\left[(\bar\eta_D - \bar\eta_d)\bar s \bar s' ++ \bar\eta_d \operatorname{diag}(\bar s)\right] . +\end{aligned} +``` + +Stare at the third line. + +The superelasticity appears in $Q$ and in $\Theta_X$, but it **cancels** from $Q + \Theta_X$. + +By {prf:ref}`mfg_prop_riccati`, the sum $Q+\Theta_X$ is the only channel through which either matrix reaches aggregate dynamics. + +Actions are the rates of change of prices, and stores pay quadratic Rotemberg costs of changing them. + +The matrix $\Gamma$ lets that cost depend on which bundle of prices is changed, with negative off-diagonal elements representing economies of scope in repricing of the kind emphasized by {cite}`Midrigan2011` and {cite}`AlvarezLippi2014`. + +There is no interaction through aggregate actions, so $\Theta_{\mathcal A} = 0$, and a nonzero rate of cost inflation makes $A$ diagonal. + +```{code-cell} ipython +def kimball(η_d, η_D, η_D_prime, s): + "Curvature and state-interaction matrices under Kimball demand." + s = np.asarray(s, dtype=float) + S = np.outer(s, s) + Q = (η_D - 1)*((η_D - η_d + η_D_prime/(η_D - 1))*S + η_d*np.diag(s)) + Θ_X = -η_D_prime*S + return Q, Θ_X + +s_bar = np.array([0.5, 0.3, 0.2]) # three products, unequal shares +η_d, η_D = 6.0, 4.0 # within-store and across-store elasticities +ρ_K, π_K = 0.04, 0.02 # discount rate and cost inflation + +n_K = len(s_bar) +A_K = π_K*np.eye(n_K) +B_K = np.eye(n_K) +Γ_K = 20.0*(np.eye(n_K) - 0.25*(np.ones((n_K, n_K)) - np.eye(n_K))) + +print("Rotemberg cost matrix Γ =\n", Γ_K) +print("\neigenvalues of Γ:", np.linalg.eigvalsh(Γ_K).round(4)) +``` + +### The superelasticity and aggregate dynamics + +Now sweep the superelasticity across a wide range, including a negative value, and watch what changes and what does not. + +```{code-cell} ipython +print(f"{'η_D′':>6} {'eigenvalues of Q':>28} {'eigenvalues of Q + Θ_X':>28}") +for η_Dp in (-3.0, 0.0, 3.0, 10.0): + Q_K, Θ_K = kimball(η_d, η_D, η_Dp, s_bar) + print(f"{η_Dp:>6} {str(np.linalg.eigvalsh(Q_K).round(3)):>28}" + f" {str(np.linalg.eigvalsh(Q_K + Θ_K).round(3)):>28}") +``` + +The own-curvature matrix $Q$ moves a great deal, and $\Theta_X$ moves with it, but the effective curvature does not move at all. + +Aggregate price dynamics are therefore invariant to the superelasticity. + +```{code-cell} ipython +print(f"{'η_D′':>6} {'aggregate eigenvalues':>34} {'individual eigenvalues':>34}") +for η_Dp in (-3.0, 0.0, 3.0, 10.0): + Q_K, Θ_K = kimball(η_d, η_D, η_Dp, s_bar) + P_K = mfg_riccati(A_K, B_K, Q_K + Θ_K, Γ_K, ρ_K) + β2_K = mfg_riccati(A_K, B_K, Q_K, Γ_K, ρ_K) + ev_agg = np.sort(np.linalg.eigvals(closed_loop(A_K, B_K, Γ_K, P_K)).real) + ev_ind = np.sort(np.linalg.eigvals(closed_loop(A_K, B_K, Γ_K, β2_K)).real) + print(f"{η_Dp:>6} {str(ev_agg.round(4)):>34} {str(ev_ind.round(4)):>34}") +``` + +The left block of numbers is identical in every row, and the right block is not. + +The superelasticity changes how an individual store behaves without changing how the industry behaves. + +Almost all of the movement is in a single mode, because $\Theta_X$ is rank one and points in the direction $\bar s$, which is the store's own share-weighted price index. + +The remaining modes, which involve relative prices within the store, barely move. + +### Static complementarity + +It is worth seeing how much static complementarity we are varying. + +Maximizing the static profit function gives the best response $x^{*}(X) = -Q^{-1}\Theta_X X$, and {cite}`AlvarezArgente2026` show that + +```{math} +:label: mfg_kimball_br +\frac{\partial x_i^{*}(X)}{\partial X_j} = \bar s_j \frac{\kappa}{1+\kappa}, +\qquad +\kappa \equiv \frac{\bar\eta_D'}{(\bar\eta_D-1)\bar\eta_D} . +``` + +```{code-cell} ipython +print(f"{'η_D′':>6}{'pass-through':>14} best response matrix, first row") +for η_Dp in (-3.0, 0.0, 3.0, 10.0): + Q_K, Θ_K = kimball(η_d, η_D, η_Dp, s_bar) + BR = -np.linalg.solve(Q_K, Θ_K) + κ = η_Dp/((η_D - 1)*η_D) + BR_closed = κ/(1 + κ)*np.outer(np.ones(n_K), s_bar) + assert np.allclose(BR, BR_closed) + print(f"{η_Dp:>6}{BR.sum(axis=1)[0]:>14.4f} {BR[0].round(4)}") + +print("\nΘ_X symmetric: ", np.allclose(Θ_K, Θ_K.T)) +print("best response symmetric:", np.allclose(BR, BR.T)) +``` + +Static pass-through runs from $-33\%$ to $+45\%$ across these rows, and it changes sign with the superelasticity, yet every one of these economies has exactly the same aggregate dynamics. + +An intuition that reads stronger static complementarity as more aggregate propagation is therefore unreliable. + +Notice also that the best response matrix is **not** symmetric, because shares differ across products, while $\Theta_X$ always is. + +Symmetry of $\Theta_X$ is what {prf:ref}`mfg_prop_planner` needs, and it survives even when the static game looks asymmetric. + +### Closed-form eigenvalues + +Because $A_K$ is a multiple of the identity and $B_K = I$, the aggregate eigenvalues have the same form as in the scalar case, one for each eigenvalue $\omega_i$ of $(\Gamma+\Theta_{\mathcal A})^{-1}(Q+\Theta_X)$: + +$$ +\lambda_i = \frac\rho2 - \sqrt{\left(\frac\rho2 + a\right)^2 + \omega_i} . +$$ + +```{code-cell} ipython +Q_K, Θ_K = kimball(η_d, η_D, 3.0, s_bar) +P_K = mfg_riccati(A_K, B_K, Q_K + Θ_K, Γ_K, ρ_K) +JG_K = closed_loop(A_K, B_K, Γ_K, P_K) + +ω = np.linalg.eigvals(np.linalg.solve(Γ_K, Q_K + Θ_K)).real +predicted = np.sort(ρ_K/2 - np.sqrt((ρ_K/2 + π_K)**2 + ω)) + +print("eigenvalues of the closed-loop matrix:", np.sort(np.linalg.eigvals(JG_K).real).round(6)) +print("from the closed-form formula: ", predicted.round(6)) +``` + +{ref}`mfg_ex8` asks what happens when inflation differs across products, so that $A$ is diagonal but not a multiple of the identity. + +### Aggregate and individual price paths + +The decomposition {eq}`mfg_decomposition` splits a store's prices into the industry average $X$ and its own deviation $z = x - X$. + +The first follows the equilibrium matrix $P$, the second the single-agent matrix $\bar\beta_2$. + +So the invariance we found should show up as identical paths for the industry and different paths for a store that is out of line. + +```{code-cell} ipython +from scipy.linalg import expm + +t_K = np.linspace(0, 8, 200) +shock = 0.10*np.ones(n_K) # ten percent above target, all products + +fig, axes = plt.subplots(1, 2, figsize=(11, 4.2)) +for k, η_Dp in enumerate((-3.0, 0.0, 3.0, 10.0)): + Q_K, Θ_K = kimball(η_d, η_D, η_Dp, s_bar) + P_K = mfg_riccati(A_K, B_K, Q_K + Θ_K, Γ_K, ρ_K) + β2_K = mfg_riccati(A_K, B_K, Q_K, Γ_K, ρ_K) + U_XX = closed_loop(A_K, B_K, Γ_K, P_K) + U_zz = closed_loop(A_K, B_K, Γ_K, β2_K) + agg = np.array([s_bar @ expm(U_XX*t) @ shock for t in t_K]) + dev = np.array([s_bar @ expm(U_zz*t) @ shock for t in t_K]) + label = f"$\\bar\\eta_D' = {η_Dp}$" + # decreasing line widths, so that four coincident curves remain visible + axes[0].plot(t_K, 100*agg, lw=6 - 1.5*k, label=label) + axes[1].plot(t_K, 100*dev, lw=2, label=label) + +axes[0].set_title('industry price index, $\\bar s\\,\'X(t)$') +axes[1].set_title("one store's deviation, $\\bar s\\,'z(t)$") +for ax in axes: + ax.set_xlabel('years') + ax.set_ylabel('percent above target') + ax.legend() +plt.tight_layout() +plt.show() +``` + +The four curves in the left panel lie exactly on top of one another. + +The four curves in the right panel do not: the higher the superelasticity, the faster a store closes a gap between its own prices and the industry's. + +Micro and macro price flexibility are governed by different objects, and only the micro one responds to the superelasticity. + +### The planner does care + +Doubling the interaction gives the planner the effective curvature $Q + 2\Theta_X = (Q + \Theta_X) - \bar\eta_D' \bar s\bar s'$, which *does* depend on the superelasticity. + +```{code-cell} ipython +print(f"{'η_D′':>6} {'eig(Q + 2Θ_X)':>26} {'planner eigenvalues':>32}") +for η_Dp in (-3.0, 0.0, 3.0, 10.0, 20.0): + Q_K, Θ_K = kimball(η_d, η_D, η_Dp, s_bar) + M = Q_K + 2*Θ_K + P_p = mfg_riccati(A_K, B_K, M, Γ_K, ρ_K) + ev = np.sort(np.linalg.eigvals(closed_loop(A_K, B_K, Γ_K, P_p)).real) + flag = "" if np.linalg.eigvalsh(M).min() > 0 else " <- not positive definite" + print(f"{η_Dp:>6} {str(np.linalg.eigvalsh(M).round(3)):>26}" + f" {str(ev.round(4)):>32}{flag}") +``` + +So the equilibrium's insensitivity to the superelasticity is special to the equilibrium. + +The gap between what a planner would do and what the industry does widens as the superelasticity rises, and in the last row the planner's effective curvature has stopped being positive definite and the returned matrix no longer stabilizes the system. + +{ref}`mfg_ex7` locates the threshold exactly, and gives it an interpretation in terms of static pass-through. + + +## Aggregate shocks and identification + +Now add a shock that hits everyone at once. + +Let $\mathcal J$ be a compensated jump process and suppose + +$$ +dx = (B\alpha - Ax)dt + \Sigma^{1/2}dW + \Upsilon \, d\mathcal J . +$$ + +With common noise the value function must keep track of the aggregate state as well as the individual state, so it becomes $v(x, X)$. + +Nevertheless, matching coefficients in the recursive HJB equation delivers a result worth emphasizing. + +```{prf:proposition} +:label: mfg_prop_common_noise + +With common noise, the matrix $P$ governing aggregate dynamics solves the **same** equilibrium Riccati equation as before. +``` + +This is a certainty-equivalence result of the kind familiar from linear quadratic control. + +It also yields a clean decomposition. + +Writing $z = x - X$ for an agent's deviation from the cross-sectional mean, + +```{math} +:label: mfg_decomposition +\begin{aligned} +dX &= \left(\Lambda P - A\right) X dt + \Upsilon d\mathcal J , \\ +dz &= \left(B\Gamma^{-1}B'\bar\beta_2 - A\right) z \, dt + \Sigma^{1/2}dW , \\ +\alpha &= (\Gamma+\Theta_{\mathcal A})^{-1}B'P X + \Gamma^{-1}B'\bar\beta_2 z . +\end{aligned} +``` + +Look carefully at the second line. + +The dynamics of an agent's deviation from the mean involve $\bar\beta_2$, the solution of the **single-agent** Riccati equation, and neither interaction matrix appears. + +The same is true of the part of the action that responds to $z$. + +```{prf:proposition} +:label: mfg_prop_identification + +Data on $x(t) - X(t)$ and $\alpha(t)-\mathcal A(t)$ contain no information about $\Theta_X$ or $\Theta_{\mathcal A}$. +``` + +This is a sharp statement of the "missing intercept" problem. + +Removing time effects from micro data -- a standard way to control for aggregate conditions -- removes exactly the variation that identifies strategic interaction. + +Aggregate data, by contrast, do carry that information, and can be combined with micro data to recover it. + +```{code-cell} ipython +# reduced-form matrices that an econometrician could estimate +U_XX = closed_loop(A, B, Γ + Θ_A, P) # drift of the aggregate state +U_zz = closed_loop(A, B, Γ, β2) # drift of the deviation from the mean + +print("aggregate drift U_XX =\n", U_XX.round(4)) +print("\ndeviation drift U_zz =\n", U_zz.round(4)) + +# now double the state interaction and recompute +P_alt = mfg_riccati(A, B, Q + 2*Θ_X, Γ + Θ_A, ρ) +print("\nwith Θ_X doubled:") +print(" aggregate drift changes: ", + not np.allclose(U_XX, closed_loop(A, B, Γ + Θ_A, P_alt))) +print(" deviation drift changes: ", + not np.allclose(U_zz, closed_loop(A, B, Γ, β2))) +``` + +{ref}`mfg_ex4` shows how to recover $\Theta_X$ from aggregate data, and why $\Theta_X$ and $\Theta_{\mathcal A}$ cannot be told apart without a normalization. + +## Exercises + +```{exercise} +:label: mfg_ex1 + +Formulas {eq}`mfg_scalar_existence` and {eq}`mfg_scalar_stability` divide the scalar model into four regions. + +Define + +$$ +-\theta^{*} \equiv q + a(a+\rho)\frac{\gamma+\theta_{\mathcal A}}{b^2}, +\qquad +-\theta^{**} \equiv q + \left(\frac\rho2+a\right)^2\frac{\gamma+\theta_{\mathcal A}}{b^2} . +$$ + +Take $\rho = 0.5$, $a = 0.3$, $b = 1$, $q = 1$, $\gamma = 1$, $\theta_{\mathcal A}=0$. + +1. Compute $\theta^{*}$ and $\theta^{**}$, and verify that $-\theta^{**} = -\theta^{*} + (\rho/2)^2(\gamma+\theta_{\mathcal A})/b^2$. +1. For values of $\theta_X$ on either side of each threshold, compute $\lambda$ and classify the outcome: stable, a unit root, divergent but admissible, or no equilibrium. +1. What does `mfg_riccati` do when no equilibrium exists? +1. Plot $\lambda$ against $\theta_X$ and mark the thresholds. +``` + +```{solution-start} mfg_ex1 +:class: dropdown +``` + +```{code-cell} ipython +ρ_e, a_e, b_e, q_e, γ_e, θA_e = 0.5, 0.3, 1.0, 1.0, 1.0, 0.0 + +θ_star = -(q_e + a_e*(a_e + ρ_e)*(γ_e + θA_e)/b_e**2) +θ_ss = -(q_e + (ρ_e/2 + a_e)**2*(γ_e + θA_e)/b_e**2) + +print(f"θ* = {θ_star:.4f}") +print(f"θ** = {θ_ss:.4f}") +print(f"gap = {θ_star - θ_ss:.4f}, " + f"(ρ/2)²(γ+θ_A)/b² = {(ρ_e/2)**2*(γ_e + θA_e)/b_e**2:.4f}") +``` + +```{code-cell} ipython +for θ_X in (0.5, -0.5, θ_star, -1.27, -1.31): + inside = (ρ_e/2 + a_e)**2 + b_e**2*(q_e + θ_X)/(γ_e + θA_e) + if inside < 0: + print(f"θ_X = {θ_X:+.4f}: no equilibrium") + continue + λ = ρ_e/2 - np.sqrt(inside) + if λ < -1e-9: + kind = "stable, X converges to zero" + elif abs(λ) <= 1e-9: + kind = "unit root, X stays where it starts" + else: + kind = "divergent but admissible" + print(f"θ_X = {θ_X:+.4f}: λ = {λ:+.5f} {kind}") +``` + +Below $\theta^{**}$ there is no real root, and the solver reports failure rather than returning a spurious answer. + +```{code-cell} ipython +try: + mfg_riccati(np.array([[a_e]]), np.array([[b_e]]), + np.array([[q_e - 1.31]]), np.array([[γ_e]]), ρ_e) +except Exception as e: + print(f"solver raises {type(e).__name__} when no equilibrium exists") +``` + +```{code-cell} ipython +θ_grid = np.linspace(-1.30, 0.5, 400) +inside = (ρ_e/2 + a_e)**2 + b_e**2*(q_e + θ_grid)/(γ_e + θA_e) +λ_grid = np.where(inside >= 0, ρ_e/2 - np.sqrt(np.maximum(inside, 0)), np.nan) + +fig, ax = plt.subplots(figsize=(8, 4.5)) +ax.plot(θ_grid, λ_grid, lw=2) +ax.axhline(0, color='k', lw=0.8) +ax.axhline(ρ_e/2, color='grey', ls=':', lw=1, label=r'$\rho/2$') +ax.axvline(θ_star, color='r', ls='--', lw=1, label=r'$\theta^{*}$') +ax.axvline(θ_ss, color='b', ls='--', lw=1, label=r'$\theta^{**}$') +ax.set_xlabel(r'$\theta_X$') +ax.set_ylabel(r'$\lambda$') +ax.legend() +plt.tight_layout() +plt.show() +``` + +Moving left along the horizontal axis means stronger complementarity in states. + +Persistence rises smoothly until $\theta^{*}$, where the aggregate state stops converging; between $\theta^{*}$ and $\theta^{**}$ equilibrium paths diverge, though slowly enough to keep discounted payoffs finite; below $\theta^{**}$ no equilibrium exists at all. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex2 + +The lecture showed that stronger complementarity in states raises the **trace** of the closed-loop matrix. + +Does it also raise every eigenvalue? + +Draw many random two-dimensional economies: take $B = I$, let $Q$ and $\Gamma$ be random positive definite matrices, let $\Theta_X$ be random negative definite, and set $\rho = 0.05$. + +For each draw, compare the equilibrium with interaction $\Theta_X$ to the one with $1.5\,\Theta_X$, keeping only draws in which both equilibria exist and are admissible. + +Report how often the trace rises and how often *every* eigenvalue rises. +``` + +```{solution-start} mfg_ex2 +:class: dropdown +``` + +```{code-cell} ipython +rng = np.random.default_rng(3) +trace_rose = eig_rose = kept = 0 + +for _ in range(400): + A_r = rng.normal(size=(2, 2))*0.3 + 0.4*np.eye(2) + B_r = np.eye(2) + M = rng.normal(size=(2, 2)); Q_r = M @ M.T + 0.5*np.eye(2) + M = rng.normal(size=(2, 2)); Γ_r = M @ M.T + 0.5*np.eye(2) + M = rng.normal(size=(2, 2)); Θ_r = -0.2 * (M @ M.T) + + try: + P0 = mfg_riccati(A_r, B_r, Q_r + Θ_r, Γ_r, 0.05) + P1 = mfg_riccati(A_r, B_r, Q_r + 1.5*Θ_r, Γ_r, 0.05) + except Exception: + continue + + J0 = closed_loop(A_r, B_r, Γ_r, P0) + J1 = closed_loop(A_r, B_r, Γ_r, P1) + if max(np.linalg.eigvals(J0).real.max(), np.linalg.eigvals(J1).real.max()) >= 0.025: + continue + + kept += 1 + trace_rose += np.trace(J1) > np.trace(J0) - 1e-10 + e0 = np.sort(np.linalg.eigvals(J0).real) + e1 = np.sort(np.linalg.eigvals(J1).real) + eig_rose += np.all(e1 >= e0 - 1e-10) + +print(f"admissible draws: {kept}") +print(f"trace rose: {trace_rose}/{kept}") +print(f"every eigenvalue rose: {eig_rose}/{kept}") +``` + +The trace result holds in every draw, as the theory says it must. + +Eigenvalue-by-eigenvalue monotonicity fails in a noticeable minority of cases. + +The reason is that $Q+\Theta_X$ and $B(\Gamma+\Theta_{\mathcal A})^{-1}B'$ need not commute, so the problem does not separate into independent one-dimensional problems, and a change that slows adjustment overall can still speed it up along some direction. + +Additional restrictions -- for instance a scalar drift matrix $A = a I$ -- restore monotonicity mode by mode. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex3 + +Compare the planner with the decentralized equilibrium. + +Using the two-dimensional example from the lecture, compute the eigenvalues of the closed-loop matrix for the equilibrium and for the planner in four cases: + +1. only complementarity in states, $\Theta_X \prec 0$ and $\Theta_{\mathcal A}=0$ +1. only substitutability in actions, $\Theta_X = 0$ and $\Theta_{\mathcal A} = 0.2 I$ +1. only complementarity in actions, $\Theta_X = 0$ and $\Theta_{\mathcal A} = -0.2 I$ +1. complementarity in both + +In which cases is the planner's allocation more persistent than the equilibrium? Explain. +``` + +```{solution-start} mfg_ex3 +:class: dropdown +``` + +```{code-cell} ipython +def eigen_pair(Θ_X_use, Θ_A_use): + "Closed-loop eigenvalues for the equilibrium and for the planner." + P_eq = mfg_riccati(A, B, Q + Θ_X_use, Γ + Θ_A_use, ρ) + P_pl = mfg_riccati(A, B, Q + 2*Θ_X_use, Γ + 2*Θ_A_use, ρ) + e_eq = np.sort(np.linalg.eigvals(closed_loop(A, B, Γ + Θ_A_use, P_eq)).real) + e_pl = np.sort(np.linalg.eigvals(closed_loop(A, B, Γ + 2*Θ_A_use, P_pl)).real) + return e_eq, e_pl + +zero = np.zeros((2, 2)) +cases = {'states, complements': (Θ_X, zero), + 'actions, substitutes': (zero, 0.2*np.eye(2)), + 'actions, complements': (zero, -0.2*np.eye(2)), + 'both complements': (Θ_X, -0.2*np.eye(2))} + +print(f"{'case':24}{'equilibrium':>22}{'planner':>22} planner slower?") +for label, (tx, ta) in cases.items(): + e_eq, e_pl = eigen_pair(tx, ta) + slower = bool(np.all(e_pl >= e_eq - 1e-12)) + print(f"{label:24}{str(e_eq.round(4)):>22}{str(e_pl.round(4)):>22} {slower}") +``` + +With only complementarity in states, the planner's allocation is more persistent. + +The planner recognizes that when one agent stays away from the steady state, others are content to stay away too, so there is less reason to hurry back. + +With only complementarity in **actions**, the comparison reverses and the planner adjusts faster: the planner internalizes that when one agent adjusts, others want to adjust as well. + +The case of substitutability in actions is the mirror image of the last one, and again makes the planner slower. + +When complementarity is present in both states and actions, the two forces work against each other. + +In this calibration the state interaction dominates and the planner is still slower, but that ranking is a quantitative accident, not a theorem: in general the comparison between planner and equilibrium persistence cannot be signed. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex4 + +This exercise works through {prf:ref}`mfg_prop_identification` and its consequences. + +1. Verify that the drift of an agent's deviation from the cross-sectional mean is unchanged when $\Theta_X$ and $\Theta_{\mathcal A}$ change. +1. Suppose an econometrician knows $\rho$, $Q$, $A$, $B$, $\Gamma$ and $\Theta_{\mathcal A}$, and estimates the aggregate drift matrix $U_{\dot X, X}$. Show how to recover $\Theta_X$, and verify your procedure numerically. +1. In the scalar model, show that $\theta_X$ and $\theta_{\mathcal A}$ are not separately identified: construct a family of pairs that generate identical aggregate dynamics *and* identical policy coefficients. +``` + +```{solution-start} mfg_ex4 +:class: dropdown +``` + +For part 1, the deviation drift is $B\Gamma^{-1}B'\bar\beta_2 - A$, and $\bar\beta_2$ solves the single-agent Riccati equation {eq}`mfg_beta2`, in which no interaction matrix appears. + +```{code-cell} ipython +for label, (tx, ta) in {'baseline': (Θ_X, Θ_A), + 'very different': (3*Θ_X, -0.1*np.eye(2))}.items(): + β2_case = mfg_riccati(A, B, Q, Γ, ρ) # does not depend on tx, ta + print(f"{label:16}: U_zz =", closed_loop(A, B, Γ, β2_case).round(6).tolist()) +``` + +For part 2, invert the definition of the closed-loop matrix to get $P$, then read $\Theta_X$ off the equilibrium Riccati equation: + +$$ +P = \Lambda^{-1}\left(U_{\dot X, X} + A\right), +\qquad +\Theta_X = P\Lambda P - Q - \rho P - PA - A'P . +$$ + +```{code-cell} ipython +U_obs = closed_loop(A, B, Γ + Θ_A, P) # what the econometrician estimates + +P_hat = np.linalg.solve(Λ, U_obs + A) +Θ_X_hat = P_hat @ Λ @ P_hat - Q - ρ*P_hat - P_hat @ A - A.T @ P_hat + +print("true Θ_X =\n", Θ_X.round(6)) +print("\nrecovered Θ_X =\n", Θ_X_hat.round(6)) +print(f"\nmaximum error: {np.abs(Θ_X - Θ_X_hat).max():.2e}") +``` + +For part 3, formula {eq}`mfg_scalar_lambda` shows that $\lambda$ depends on the two interactions only through the ratio $(q+\theta_X)/(\gamma+\theta_{\mathcal A})$. + +Holding that ratio fixed traces out a one-parameter family of observationally equivalent economies. + +The policy coefficient cannot break the tie, because $\lambda = b\,U_{\alpha,X} - a$ ties it to $\lambda$. + +```{code-cell} ipython +ratio = (q + (-0.4))/(γ + 0.3) # baseline θ_X = -0.4, θ_A = 0.3 + +print(f"{'θ_A':>7}{'θ_X':>12}{'λ':>12}{'policy coeff':>15}") +for θ_A_alt in (0.3, 0.0, 0.6, 1.0): + θ_X_alt = ratio*(γ + θ_A_alt) - q + λ_alt = scalar_lambda(ρ_s, a, b, q, γ, θ_X_alt, θ_A_alt) + print(f"{θ_A_alt:>7.2f}{θ_X_alt:>12.6f}{λ_alt:>12.6f}{(λ_alt + a)/b:>15.6f}") +``` + +Every row describes a different economy, with a different amount of strategic complementarity in states and in actions, yet all of them generate exactly the same aggregate dynamics and the same policy rule. + +Distinguishing them requires normalizing one interaction, or bringing outside information to bear. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex5 + +This exercise derives the two curvatures that we simply assumed when calibrating the capital accumulation example. + +Suppose the adjustment cost is + +$$ +\psi(i) = \frac{\phi}{2}\,\bar i\left(\frac{i - \bar i}{\bar i}\right)^2 , +$$ + +so that $\psi'(\bar i) = 0$ and $\psi''(\bar i) = \phi/\bar i$, and suppose investment goods are supplied with constant elasticity $\varepsilon_s$, + +$$ +\mathcal P(I) = \left(\frac{I}{\bar i}\right)^{1/\varepsilon_s} . +$$ + +1. Use the steady-state Euler equation $\Pi_k(\bar k, \bar k) = \rho + \delta$ to show that $\bar k^{1-\nu} = \nu(\eta-1)/[\eta(\rho+\delta)]$, and then that + +$$ +\gamma = \phi\,\delta\bar k^{1-\nu} , \qquad +\theta_{\mathcal A} = \frac{\delta \bar k^{1-\nu}}{\varepsilon_s} . +$$ + +1. Which $(\phi, \varepsilon_s)$ reproduce the calibration $\gamma = 0.05$ and $\theta_{\mathcal A} = 0.09$? +1. Show that all pairs with the same value of $\phi + 1/\varepsilon_s$ generate identical *aggregate* dynamics, and verify numerically that they nevertheless imply different individual investment rules and different planner allocations. +``` + +```{solution-start} mfg_ex5 +:class: dropdown +``` + +For part 1, $\Pi_k(k,K) = \nu\frac{\eta-1}{\eta}k^{\nu(1-1/\eta)-1}K^{\nu/\eta}$, so on the diagonal $\Pi_k(\bar k,\bar k) = \nu\frac{\eta-1}{\eta}\bar k^{\nu-1}$. + +Setting this equal to $\rho+\delta$, which is the steady-state user cost when $\mathcal P(\bar i) = 1$ and $\psi'(\bar i) = 0$, gives the stated expression for $\bar k^{1-\nu}$. + +Since $\bar\Pi = \bar k^{\nu}$ and $\bar i = \delta\bar k$, + +$$ +\gamma = \frac{\delta^2\bar k^2}{\bar k^{\nu}}\frac{\phi}{\delta \bar k} += \phi\,\delta\bar k^{1-\nu} , +\qquad +\theta_{\mathcal A} = \frac{\delta^2\bar k^2}{\bar k^{\nu}} +\frac{1}{\varepsilon_s \delta\bar k} = \frac{\delta\bar k^{1-\nu}}{\varepsilon_s} . +$$ + +```{code-cell} ipython +scale_c = δ_c*ν_c*(η_c - 1)/(η_c*(ρ_c + δ_c)) # δ k̄^{1-ν} + +print(f"δ k̄^(1-ν) = {scale_c:.4f}") +print(f"φ implied by γ = {γ_c}: {γ_c/scale_c:.4f}") +print(f"ε_s implied by θ_A = {θ_A_c}: {scale_c/θ_A_c:.4f}") +``` + +The calibration corresponds to a moderately convex adjustment cost, $\phi = 0.14$, and an investment supply elasticity of about $3.9$. + +For part 3, $\gamma + \theta_{\mathcal A} = \delta\bar k^{1-\nu}(\phi + 1/\varepsilon_s)$, and by {eq}`mfg_scalar_lambda` only this sum matters for $\lambda$. + +But the individual's Riccati equation {eq}`mfg_beta2` involves $\gamma$ alone, and the planner's involves $\gamma + 2\theta_{\mathcal A}$, so both separate the pairs. + +```{code-cell} ipython +print(f"{'φ':>6}{'1/ε_s':>8}{'aggregate':>12}{'individual':>12}{'planner':>10}") +for φ, inv_ε in ((0.40, 0.00), (0.30, 0.10), (1/7, 0.257143), (0.05, 0.35)): + γ_case, θ_A_case = φ*scale_c, inv_ε*scale_c + λ_agg = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_case, θ_X_c, θ_A_case) + λ_ind = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_case, 0.0, 0.0) + λ_pl = scalar_lambda(ρ_c, δ_c, δ_c, q_c, γ_case, 2*θ_X_c, 2*θ_A_case) + print(f"{φ:>6.3f}{inv_ε:>8.3f}{λ_agg:>12.6f}{λ_ind:>12.6f}{λ_pl:>10.4f}") +``` + +Aggregate dynamics are identical down to the last digit, while individual behavior ranges from sluggish to almost frictionless. + +An econometrician with only aggregate data cannot tell a congested market for capital goods from an expensive installation technology, which is the scalar version of {prf:ref}`mfg_prop_identification`. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex6 + +The right panel of the figure above suggests that something special happens as returns to scale approach constancy. + +1. Show that $q + \theta_X \to 0$ and $\lambda \to -\delta$ as $\nu \to 1$, for every $\eta > 1$. +1. Show that the response of aggregate investment to a capital gap, $(\lambda+\delta)/\delta$, goes to zero in the same limit. +1. Explain the limit in economic terms. +``` + +```{solution-start} mfg_ex6 +:class: dropdown +``` + +Part 1 is immediate from {eq}`mfg_capital_coeffs`: $q + \theta_X = \frac{\eta-1}{\eta}\nu(1-\nu) \to 0$, and with $a = b = \delta$, + +$$ +\lambda \to \frac\rho2 - \sqrt{\left(\frac\rho2 + \delta\right)^2} = -\delta . +$$ + +For part 2, the scalar version of {eq}`mfg_aggregate_action` gives $\mathcal A = \frac{\lambda+\delta}{\delta}X$, which vanishes as $\lambda \to -\delta$. + +```{code-cell} ipython +print(f"{'ν':>8}{'q + θ_X':>10}{'λ':>10}{'half-life':>11}{'A/X':>9}") +for ν in (0.7, 0.9, 0.99, 0.999, 0.99999): + λ_ν = scalar_lambda(ρ_c, δ_c, δ_c, cap_q(η_c, ν), γ_c, cap_θ_X(η_c, ν), θ_A_c) + print(f"{ν:>8}{cap_q(η_c,ν) + cap_θ_X(η_c,ν):>10.5f}{λ_ν:>10.5f}" + f"{half_life(λ_ν):>11.3f}{(λ_ν + δ_c)/δ_c:>9.4f}") +``` + +For part 3, note that $q + \theta_X = -\bar k^2\left[\Pi_{kk} + \Pi_{kK}\right]/\bar\Pi$ measures the curvature of profit along the *diagonal*, that is, the rate at which the marginal profitability of capital falls when the whole industry expands together. + +With $\nu = 1$ the profit function {eq}`mfg_capital_profit` is homogeneous of degree one in $(k, K)$ jointly, so $\Pi_k(k,k)$ does not depend on $k$ at all. + +An industry-wide capital gap then creates no incentive to invest differently, and the gap closes only through depreciation. + +Effective curvature is zero, the economy sits exactly on the boundary of the existence region in {eq}`mfg_scalar_existence`, and the equilibrium remains unique and stable with $\lambda = -\delta$. + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex7 + +In the Kimball example, the equilibrium's effective curvature $Q + \Theta_X$ is positive definite for every superelasticity, but the planner's, $Q + 2\Theta_X$, is not. + +1. Using $Q + 2\Theta_X = (Q+\Theta_X) - \bar\eta_D'\,\bar s\bar s'$, show that it is positive definite if and only if +$\bar\eta_D'\,\bar s'(Q+\Theta_X)^{-1}\bar s < 1$. +1. Show that $(Q+\Theta_X)\mathbf 1 = (\bar\eta_D-1)\bar\eta_D\,\bar s$, where $\mathbf 1$ is a vector of ones, and hence that the condition is $\bar\eta_D' < (\bar\eta_D-1)\bar\eta_D$. +1. Show that this is exactly the condition that static pass-through in {eq}`mfg_kimball_br` be below one half, and verify the threshold numerically. +``` + +```{solution-start} mfg_ex7 +:class: dropdown +``` + +For part 1, if $\bar\eta_D' \leq 0$ the rank-one term is added rather than subtracted and positive definiteness is immediate. + +If $\bar\eta_D' > 0$, then for $M \succ 0$ the matrix $M - c vv'$ with $c>0$ is positive definite if and only if $c\, v'M^{-1}v < 1$, which follows from the determinant identity $\det(M - cvv') = \det(M)(1 - c\,v'M^{-1}v)$ applied to every leading block, or directly from the Schur complement of the bordered matrix. + +For part 2, using $\bar s'\mathbf 1 = 1$ in the third line of {eq}`mfg_kimball_Q`, + +$$ +(Q+\Theta_X)\mathbf 1 += (\bar\eta_D-1)\left[(\bar\eta_D-\bar\eta_d)\bar s + \bar\eta_d \bar s\right] += (\bar\eta_D-1)\bar\eta_D\,\bar s . +$$ + +So $(Q+\Theta_X)^{-1}\bar s = \mathbf 1/[(\bar\eta_D-1)\bar\eta_D]$ and therefore + +$$ +\bar s'(Q+\Theta_X)^{-1}\bar s = \frac{1}{(\bar\eta_D-1)\bar\eta_D} , +$$ + +which turns the condition in part 1 into $\bar\eta_D' < (\bar\eta_D-1)\bar\eta_D$. + +For part 3, that inequality says precisely that $\kappa < 1$, and $\kappa/(1+\kappa)$ is increasing in $\kappa$ with value $1/2$ at $\kappa = 1$. + +Total static pass-through is $\sum_j \partial x_i^{*}/\partial X_j = \kappa/(1+\kappa)$ because the shares sum to one, so the planner's problem is well behaved exactly when a store would pass less than half of an industry-wide price increase into its own prices. + +```{code-cell} ipython +lo, hi = 0.0, 100.0 +for _ in range(60): + mid = (lo + hi)/2 + Q_m, Θ_m = kimball(η_d, η_D, mid, s_bar) + if np.linalg.eigvalsh(Q_m + 2*Θ_m).min() > 0: + lo = mid + else: + hi = mid + +κ_star = lo/((η_D - 1)*η_D) +print(f"threshold by bisection: η_D′ = {lo:.6f}") +print(f"(η_D - 1) η_D = {(η_D - 1)*η_D:.6f}") +print(f"κ at the threshold = {κ_star:.6f}") +print(f"pass-through there = {κ_star/(1 + κ_star):.6f}") +``` + +```{solution-end} +``` + +```{exercise} +:label: mfg_ex8 + +The closed-form eigenvalues in the Kimball example used $A = \pi I$. + +Replace it by $A = \operatorname{diag}(0.00, 0.02, 0.06)$, so that the three products face different rates of cost inflation. + +1. Check that the closed-form formula fails. +1. Check that the trace comparative static of the **Persistence** section still holds, by scaling $\Theta_X$ and recording the trace and the eigenvalues of the closed-loop matrix. +1. Does the invariance of aggregate dynamics to the superelasticity survive? +``` + +```{solution-start} mfg_ex8 +:class: dropdown +``` + +```{code-cell} ipython +A_het = np.diag([0.00, 0.02, 0.06]) +Q_h, Θ_h = kimball(η_d, η_D, 3.0, s_bar) + +P_h = mfg_riccati(A_het, B_K, Q_h + Θ_h, Γ_K, ρ_K) +JG_h = closed_loop(A_het, B_K, Γ_K, P_h) + +ω_h = np.linalg.eigvals(np.linalg.solve(Γ_K, Q_h + Θ_h)).real +print("closed-loop eigenvalues:", np.sort(np.linalg.eigvals(JG_h).real).round(5)) +print("closed-form formula: ", + np.sort(ρ_K/2 - np.sqrt((ρ_K/2 + np.diag(A_het))**2 + np.sort(ω_h))).round(5)) +``` + +The formula no longer holds: with $A$ not a multiple of the identity there is no single scalar shift to put inside the square root, and the eigenvectors of $A$ and of $Q+\Theta_X$ need not agree. + +The errors are modest here because the three inflation rates are close together, but they are systematic, and they grow with the dispersion in $A$. + +The trace result, on the other hand, requires no such restriction. + +```{code-cell} ipython +print(f"{'scale on Θ_X':>14}{'trace':>10} eigenvalues") +for scale in (0.0, 0.5, 1.0, 1.5): + P_s = mfg_riccati(A_het, B_K, Q_h + scale*Θ_h, Γ_K, ρ_K) + JG_s = closed_loop(A_het, B_K, Γ_K, P_s) + print(f"{scale:>14}{np.trace(JG_s):>10.5f} " + f"{np.sort(np.linalg.eigvals(JG_s).real).round(5)}") +``` + +Raising the scale strengthens complementarity and raises the trace, exactly as in the two-dimensional example earlier. + +For part 3, the invariance has nothing to do with $A$: it comes from the cancellation of $\bar\eta_D'$ in $Q + \Theta_X$, and {prf:ref}`mfg_prop_riccati` shows that only $Q+\Theta_X$ and $\Gamma+\Theta_{\mathcal A}$ enter the equilibrium Riccati equation. + +```{code-cell} ipython +for η_Dp in (-3.0, 0.0, 3.0, 10.0): + Q_i, Θ_i = kimball(η_d, η_D, η_Dp, s_bar) + P_i = mfg_riccati(A_het, B_K, Q_i + Θ_i, Γ_K, ρ_K) + ev = np.sort(np.linalg.eigvals(closed_loop(A_het, B_K, Γ_K, P_i)).real) + print(f"η_D′ = {η_Dp:>5}: {ev.round(6)}") +``` + +```{solution-end} +``` + + +## Further reading + +{cite}`AlvarezArgente2026` develop the results in this lecture, along with the two economic examples that we implemented above. + +They also extend the analysis to interactions through higher moments of the cross-sectional distribution. + +{cite}`CarmonaDelarue2018` give a comprehensive probabilistic treatment of mean field games, and {cite}`AchdouEtAl2022` describe the numerical methods used to solve the coupled partial differential equations when the model is not linear quadratic. diff --git a/lectures/lucas_prescott_investment.md b/lectures/lucas_prescott_investment.md new file mode 100644 index 000000000..cc406f85e --- /dev/null +++ b/lectures/lucas_prescott_investment.md @@ -0,0 +1,989 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: Python 3 + language: python + name: python3 +--- + +(lucas_prescott_investment)= +```{raw} jupyter + +``` + +# Investment Under Uncertainty + +```{contents} Contents +:depth: 2 +``` + +In addition to what's in Anaconda, this lecture will need the following libraries: + +```{code-cell} ipython +--- +tags: [hide-output] +--- +!pip install quantecon +``` + +## Overview + +This lecture studies {cite}`Lucas_Prescott_1971`, a paper that helped to ignite a *rational expectations revolution*. + +Lucas and Prescott studied a competitive industry in which + +* demand shifts randomly each period +* firms face costs of adjusting their capital stocks +* firms must forecast future prices in order to decide how much to invest +* the probability distribution that firms use to forecast prices **equals** the probability distribution that their investment decisions actually generate + +That last bullet point is what {cite}`muth1961` called **rational expectations**. + +The QuantEcon lecture {doc}`rational_expectations` presents what we might call a "baby" version of the Lucas-Prescott model. + +That lecture studies a linear-quadratic industry **without uncertainty**. + +The present lecture describes the more ambitious structure that Lucas and Prescott actually built. + +Relative to the baby version, Lucas and Prescott + +* let demand be shifted by a Markov process $\{u_t\}$, so that the equilibrium is a stochastic process rather than a deterministic path +* allow a nonlinear technology for converting investment into capacity +* prove that a competitive equilibrium **exists** and is **unique** +* show that the equilibrium is the solution of a **planning problem** that maximizes discounted consumer surplus +* show that the equilibrium is a **Markov process** in the state $(k_t, u_t)$ +* provide conditions under which that Markov process has an **invariant probability distribution** to which it converges from any initial condition + +The last item is the part of the paper that most influenced later work. + +A model whose equilibrium is a Markov process with an invariant distribution is a model that can be taken to time series data. + +That observation set the stage for the *rational expectations econometrics* subsequently developed by Lars Peter Hansen and Thomas Sargent {cite}`HansenSargent1980`. + +Along the way we'll describe how {cite}`PrescottMehra1980` later distilled the Lucas-Prescott structure into a general definition of a **recursive competitive equilibrium**. + +The sequel {doc}`optimal_growth_uncertainty` studies a paper that asks the same questions about a one-sector optimal growth model, {cite}`BrockMirman1972`, and a paper that uses that model to think about Tobin's $q$, {cite}`Sargent1980q`. + +Let's start with some imports: + +```{code-cell} ipython +import numpy as np +import matplotlib.pyplot as plt +import quantecon as qe +from collections import namedtuple +``` + +## The industry + +An industry consists of many small firms. + +Each firm produces a single output $q_t$ with a single input, capital $k_t$. + +Production has constant returns to scale, and with a suitable choice of units the production function is + +```{math} +:label: lp_production +0 \leq q_t \leq k_t . +``` + +Because capital is the only input and output is sold at a positive price, every firm produces at capacity, so that $q_t = k_t$. + +Let $x_t$ denote gross investment. + +Capacity next period is related to capacity this period and investment by + +```{math} +:label: lp_accumulation +k_{t+1} = k_t \, h\!\left(\frac{x_t}{k_t}\right), +``` + +where $h$ is bounded, continuously differentiable, increasing, and strictly **concave**. + +The strict concavity of $h$ is what creates **costs of adjustment**: doubling the rate of investment per unit of capital less than doubles the resulting increment to capacity. + +Adjustment costs are why firms change their capital stocks gradually instead of jumping immediately to a long-run target. + +Assume that $\delta = h^{-1}(1)$ exists and satisfies $0 < \delta < 1$. + +Then $x_t = \delta k_t$ is the investment rate that just maintains capacity, so $\delta$ plays the role of a depreciation rate. + +Let $p_t$ be the output price and let $\beta = 1/(1+r)$, where $r > 0$ is the cost of capital. + +The present value of a firm is + +```{math} +:label: lp_value +V = \sum_{t=0}^\infty \beta^t \left[ p_t q_t - x_t \right]. +``` + +Because the allocation of a given industry capital stock across firms does not matter, we can use $k_t, x_t, q_t$ interchangeably for firm and industry variables. + +Equivalently, we can think of a competitive industry with a single price-taking firm. + +Industry demand is subject to random shifts: + +```{math} +:label: lp_demand +p_t = D(q_t, u_t), +``` + +where $D$ is continuous and strictly decreasing in $q_t$ and increasing in $u_t$, so that an increase in $u_t$ shifts the demand curve to the right. + +The demand shifter $\{u_t\}$ is a Markov process with transition distribution $p(\cdot, u)$, meaning that the probability that $u_{t+1} \in A$ conditional on $u_t = u$ is $\int_A p(dz, u)$. + +## The firm and the price of installed capital + +Before studying equilibrium, Lucas and Prescott pause over an instructive question: what does an individual firm actually need to know? + +Let $w_t$ be the current market value of a unit of installed capital and let $w^*_t$ be the value per unit expected to prevail next period. + +A firm that begins period $t$ with capital $k_t$ and invests $x$ obtains next-period capital worth $\beta k_t h(x/k_t) w^*_t$ at a cost of $x$, so it solves + +$$ +\max_x \left[ -x + \beta k_t h(x/k_t) w^*_t \right] . +$$ + +The current value of the firm is + +```{math} +:label: lp_firmvalue +w_t k_t = p_t k_t - x + \beta k_t h(x/k_t) w^*_t +``` + +and the first-order condition is + +```{math} +:label: lp_firmfoc +0 \geq -1 + \beta h'(x/k_t) w^*_t, \quad \text{with equality if } x > 0 . +``` + +Solving {eq}`lp_firmvalue` and {eq}`lp_firmfoc` jointly for $x$ and $w^*_t$ gives an investment function of the form + +```{math} +:label: lp_investment_fn +x_t = k_t \, g(w_t - p_t), \qquad g'(\cdot) > 0 . +``` + +Lucas and Prescott note that {eq}`lp_investment_fn` is essentially the investment function that Grunfeld had used in empirical work, with the market value of the firm as an explanatory variable. + +Their argument is stronger than Grunfeld's, though: a firm does not need to forecast its own future income stream at all. + +It needs only to know the value that securities markets place on a unit of installed capital. + +Readers will recognize a version of what later became known as Tobin's $q$ theory of investment. + +But equation {eq}`lp_investment_fn` is a **consistency requirement**, not yet a theory of capital accumulation, because the path of $w_t$ is still unknown. + +To determine $w_t$ we have to study equilibrium. + +## Rational expectations equilibrium + +Firms must forecast future prices. + +Lucas and Prescott describe the usual approach as postulating a forecasting rule -- for example "adaptive expectations" -- that generates investment behavior, which in conjunction with demand generates an actual price process. + +They object that if the underlying disturbance has a regular stochastic character, then, except by coincidence, forecast prices and actual prices will have **different probability distributions**, and the difference will be persistent, costly, and easy to correct. + +So they go to the opposite extreme and assume that the actual and anticipated prices have the **same probability distribution**. + +To say this precisely, fix an initial state $(k_0, u_0)$. + +Because prices depend on the history of demand shocks, an anticipated price process is a sequence $\{p_t\}$ of functions of $(u_1, \ldots, u_t)$. + +Similarly, an investment-output plan is a pair of sequences $\{q_t, x_t\}$ of functions of $(u_1, \ldots, u_t)$ -- a contingency plan that says in advance what the firm will do after every possible history. + +```{prf:definition} +:label: lp_equilibrium_def + +An **industry equilibrium** for a fixed initial state $(k_0, u_0)$ is a triple of sequences $\{q^0_t, x^0_t, p^0_t\}$ such that + +1. the demand curve {eq}`lp_demand` holds for every history, and +1. the plan $\{q^0_t, x^0_t\}$ maximizes expected present value + + $$ + E \left\{ \sum_{t=0}^\infty \beta^t \left[ p^0_t q_t - x_t \right] \right\} + $$ + + over all plans $\{q_t, x_t\}$ that satisfy {eq}`lp_production` and {eq}`lp_accumulation`, **given** the price process $\{p^0_t\}$. +``` + +The rational expectations requirement is hiding in plain sight in {prf:ref}`lp_equilibrium_def`. + +The price process $\{p^0_t\}$ that firms take as given when they maximize is the *same* price process that their own decisions generate through the demand curve. + +This is exactly the fixed point idea of the lecture {doc}`rational_expectations`, where a **perceived law of motion** $H$ for aggregate output must equal the **actual law of motion** that the resulting decision rule generates. + +The difference is that here the fixed point is in a space of sequences of functions of histories rather than in a space of linear decision rules. + +```{note} +Lucas and Prescott are careful about what rationality does and does not assume. + +They write that they "surrender, in advance, any hope of shedding light on the process by which firms translate current information into price forecasts." + +They also defend the assumption: if the demand shift process really does have a regular, stationary structure, then expectations that are rational in their sense "are surely more plausible than any simple, adaptive scheme"; and if it does not, then adopting some other expectations hypothesis "will certainly not improve matters." +``` + +## Equilibrium as a planning problem + +How can we compute an object defined by a fixed point in such a large space? + +Lucas and Prescott's answer is the device that the lecture {doc}`rational_expectations` also uses: find a **planning problem** whose solution is the equilibrium. + +Define **consumer surplus** as the area under the demand curve + +$$ +s(q, u) = \int_0^q D(z, u) \, dz , +$$ + +and define discounted consumer surplus net of investment costs by + +```{math} +:label: lp_surplus +S = E \left\{ \sum_{t=0}^\infty \beta^t \left[ s(q_t, u_t) - x_t \right] \right\} . +``` + +Associated with the problem of maximizing $S$ is the functional equation + +```{math} +:label: lp_bellman +v(k, u) = \sup_{x \geq 0} \left\{ s(k, u) - x + \beta \int v\!\left[ k h\!\left(\frac{x}{k}\right), z \right] p(dz, u) \right\} . +``` + +```{prf:theorem} +:label: lp_theorem1 + +The functional equation {eq}`lp_bellman` has a unique bounded solution $v$, and for each $(k,u)$ the supremum is attained at a unique $x(k,u)$. + +In terms of that policy function, the unique industry equilibrium, given $(k_0, u_0)$, is + +$$ +x_t = x(k_t, u_t), \qquad +k_{t+1} = k_t h\!\left(\frac{x(k_t,u_t)}{k_t}\right), \qquad +q_t = k_t, \qquad +p_t = D(q_t, u_t) . +$$ +``` + +The proof has two halves, and both halves matter for later work. + +The first half shows that a competitive equilibrium maximizes $S$, and conversely. + +This is an application of the welfare theorems in an infinite-dimensional commodity space, using the valuation equilibria of Debreu and the price systems that Prescott and Lucas developed in a companion paper. + +Lucas and Prescott are explicit that they use this connection only as a computational device: "the welfare significance of $S$ is not important. We are interested only in using the connection between the maximization of $S$ and competitive equilibrium in order to determine the properties of the latter." + +The second half shows that the planning problem is solved by the functional equation {eq}`lp_bellman`. + +Here they use the operator + +$$ +Tf(k,u) = \sup_{x \geq 0} \left\{ s(k,u) - x + \beta \int f\!\left[ k h(x/k), z \right] p(dz, u) \right\} +$$ + +and verify that $T$ is monotone and satisfies a discounting property, so that by Blackwell's theorem {cite}`Blackwell1965` it has a unique fixed point that successive approximations converge to. + +They also show that $T$ preserves concavity and monotonicity in $k$, which delivers a unique and continuous policy function $x(k,u)$. + +```{note} +These arguments are now standard and are treated at length in {cite}`StokeyLucas1989`. + +In 1971 they were not standard, which is one reason the paper is hard to read. + +Much of its length is devoted to measurability details -- Baire functions, Borel sets -- that a modern treatment would relegate to an appendix. +``` + +Two features of {prf:ref}`lp_theorem1` deserve emphasis. + +First, the equilibrium is **recursive**: the pair $(k_t, u_t)$ is a Markov process, and equilibrium prices and quantities are time-invariant functions of it. + +Second, the equilibrium is computed **without ever iterating on a mapping from beliefs to outcomes**. + +The lecture {doc}`rational_expectations` explains why that matters: the mapping $\Phi$ from a perceived law of motion to an actual law of motion is **not a contraction**, and iterating on it can diverge. + +The planning problem replaces an unreliable fixed point calculation with a dynamic program that is a contraction. + +## Recursive competitive equilibrium + +{cite}`PrescottMehra1980` later extracted the general structure that Lucas and Prescott had exploited. + +Their goal was to replace a search for equilibrium *sequences of contingency functions*, in the style of Arrow and Debreu, with a search for equilibrium **decision rules**. + +Such rules specify current actions as functions of a small number of **state variables** that summarize the effects of past decisions and current information. + +As Prescott and Mehra put it, these equilibrium decision rules "must be time invariant in order to apply standard time series methods and this necessitates a recursive structure." + +That sentence is the bridge from Lucas and Prescott's theory to econometrics. + +In a recursive competitive equilibrium + +* the state variables should be of minimal dimension, indexing only the factors that can change over time +* the state is observed, or is an invertible function of observables +* the conditional distribution of next period's state given current decisions and the current state is time invariant +* individual decision rules are optimal given equilibrium pricing functions, and markets clear + +Prescott and Mehra note that their structure "subsumes the structure considered in Lucas and Prescott's analysis of equilibrium investment under uncertainty." + +Their analysis also establishes optimality of recursive equilibria and supportability of Pareto optima "in a simpler and more direct way" than arguments that pass through equivalence with state-contingent equilibria. + +For us, the important point is that {prf:ref}`lp_theorem1` produces exactly the objects that rational expectations econometrics needs: time-invariant decision rules, driven by a Markov state, with cross-equation restrictions linking the parameters of the shock process to the parameters of the decision rules. + +## A computable version + +Let's now compute equilibria of a version of the model. + +We take the adjustment technology + +$$ +h(z) = (1 - \delta + z)^\alpha, \qquad 0 < \alpha \leq 1 , +$$ + +which satisfies the Lucas-Prescott assumptions: $h$ is increasing and strictly concave for $\alpha < 1$, and $h(\delta) = 1$, so $\delta$ is the maintenance investment rate. + +When $\alpha = 1$ we recover the familiar linear accumulation equation $k_{t+1} = (1-\delta) k_t + x_t$. + +When $\alpha < 1$ there are adjustment costs. + +Note that $h'(\delta) = \alpha$, a fact we will use below. + +We take a linear inverse demand curve + +$$ +D(q, u) = a_0 + u - a_1 q , +$$ + +which matches the demand curve in the lecture {doc}`rational_expectations` except that it is shifted by $u$. + +Consumer surplus is then + +$$ +s(k, u) = (a_0 + u) k - \frac{a_1}{2} k^2 . +$$ + +The demand shifter follows a Gaussian AR(1) process + +$$ +u_{t+1} = \rho u_t + \sigma \epsilon_{t+1}, \qquad \epsilon_{t+1} \sim N(0,1), +$$ + +which is an example that Lucas and Prescott themselves offer of a process satisfying their assumptions. + +We discretize it with the Tauchen method. + +Rather than choosing investment $x$ directly, it is convenient to let the planner choose next period's capital $k'$ on a grid, and to invert {eq}`lp_accumulation` to find the required investment + +$$ +x = k \left[ \left(\frac{k'}{k}\right)^{1/\alpha} - (1 - \delta) \right] . +$$ + +```{code-cell} ipython +Model = namedtuple("Model", "r β δ α a0 a1 k u P X feasible s") + +def create_model(r=0.05, δ=0.10, α=0.70, a0=1.0, a1=0.01, + ρ=0.9, σ=0.02, n_u=9, n_k=400, k_lo=20.0, k_hi=160.0): + "Discretize the Lucas-Prescott industry." + β = 1 / (1 + r) + mc = qe.markov.tauchen(n_u, ρ, σ) + u, P = mc.state_values, mc.P + k = np.linspace(k_lo, k_hi, n_k) + # investment needed to move from k (rows) to k' (columns) + X = k[:, None] * ((k[None, :] / k[:, None])**(1/α) - (1 - δ)) + feasible = X >= 0 + s = (a0 + u[None, :]) * k[:, None] - a1 * k[:, None]**2 / 2 + return Model(r, β, δ, α, a0, a1, k, u, P, X, feasible, s) +``` + +We solve the planner's Bellman equation {eq}`lp_bellman` by value function iteration, with Howard policy improvement steps to speed convergence. + +```{code-cell} ipython +def solve_model(m, tol=1e-8, maxit=1000, howard=30): + "Solve the planning problem; return value function and policies." + n_k, n_u = len(m.k), len(m.u) + R = np.where(m.feasible, -m.X, -1e12) + rows, cols = np.arange(n_k)[:, None], np.arange(n_u)[None, :] + v = m.s.copy() + + for it in range(maxit): + EV = v @ m.P.T # EV[k', u] = E[v(k', u') | u] + obj = R[:, :, None] + m.β * EV[None, :, :] # (k, k', u) + idx = obj.argmax(axis=1) # choice of k' given (k, u) + v_new = m.s + np.take_along_axis(obj, idx[:, None, :], axis=1)[:, 0, :] + + for _ in range(howard): # policy evaluation steps + EV = v_new @ m.P.T + v_new = m.s + R[rows, idx] + m.β * EV[idx, cols] + + if np.max(np.abs(v_new - v)) < tol: + v = v_new + break + v = v_new + + k_next = m.k[idx] + x = np.take_along_axis(m.X, idx, axis=1) + return v, idx, k_next, x + +m = create_model() +v, idx, k_next, x = solve_model(m) +print(f"grid: {len(m.k)} capital points, {len(m.u)} demand states") +``` + +Let's look at the equilibrium investment policy and the law of motion for capital. + +```{code-cell} ipython +fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) + +for j in [0, len(m.u)//2, len(m.u)-1]: + axes[0].plot(m.k, x[:, j], label=f'$u = {m.u[j]:.3f}$') + axes[1].plot(m.k, k_next[:, j], label=f'$u = {m.u[j]:.3f}$') + +axes[0].plot(m.k, m.δ * m.k, 'k--', lw=1, label=r'$\delta k$') +axes[0].set_xlabel('$k$'); axes[0].set_ylabel('$x(k, u)$') +axes[0].set_title('investment policy') +axes[1].plot(m.k, m.k, 'k--', lw=1, label='45 degree line') +axes[1].set_xlabel('$k$'); axes[1].set_ylabel("$k'(k, u)$") +axes[1].set_title('law of motion for capital') +for ax in axes: + ax.legend() +plt.tight_layout() +plt.show() +``` + +Capital rises when $x(k,u)$ lies above the maintenance line $\delta k$ and falls when it lies below. + +Higher demand shifts the investment policy up, so the capital stock that the industry sustains is higher when demand is strong. + +### The price of installed capital + +The planner's problem also delivers the market value $w$ of a unit of installed capital that appeared in {eq}`lp_firmvalue`. + +Let $z = x/k$ denote the investment rate. + +Differentiating the Bellman equation {eq}`lp_bellman` and using the envelope condition gives + +$$ +w(k,u) = v_k(k, u) = D(k,u) + \beta E\left[ v_k(k', u') \mid u \right] \left[ h(z) - z h'(z) \right], +$$ + +while the first-order condition for $x$ is + +```{math} +:label: lp_planner_foc +\beta E\left[ v_k(k', u') \mid u \right] = \frac{1}{h'(z)} . +``` + +Combining them expresses the shadow price in closed form, + +```{math} +:label: lp_shadow +w(k,u) = D(k,u) + \frac{h(z)}{h'(z)} - z . +``` + +The marginal value of installed capital equals the current price of output plus the value of the capacity that the unit carries into the future. + +Notice that {eq}`lp_planner_foc` is the planner's counterpart of the firm's first-order condition {eq}`lp_firmfoc`, with $w^* = E[v_k(k',u') \mid u]$. + +That correspondence is the "big $K$, little $k$" logic of the lecture {doc}`rational_expectations` in its Lucas-Prescott form. + +```{code-cell} ipython +h = lambda z, m: (1 - m.δ + z)**m.α +h_prime = lambda z, m: m.α * (1 - m.δ + z)**(m.α - 1) + +z = x / m.k[:, None] +D = m.a0 + m.u[None, :] - m.a1 * m.k[:, None] +w = D + h(z, m) / h_prime(z, m) - z + +# check the first-order condition (up to grid error) +Ew = np.take_along_axis(w @ m.P.T, idx, axis=0) +resid = np.abs(m.β * Ew - 1 / h_prime(z, m)) +scale = np.median(1 / h_prime(z, m)) +near = (m.k > 70) & (m.k < 95) # capital levels the industry actually visits +print(f"typical size of each side of the FOC: {scale:.3f}") +print(f"median residual, all k: {np.median(resid):.2e}") +print(f"median residual, 70 < k < 95: {np.median(resid[near]):.2e}") +``` + +The residual is a few tenths of one per cent of the magnitude of the terms being compared, which confirms {eq}`lp_shadow`. + +It does not vanish entirely because the planner chooses $k'$ from a finite grid, so the policy function jumps by a whole grid point at a time. + +Residuals are much larger at capital levels so extreme that the industry never visits them. + +## Long run behavior with serially independent demand + +Lucas and Prescott next ask what happens in the long run. + +They treat two cases, and the first is the special case in which $u_t$ and $u_s$ are independent for $s \neq t$. + +Inspecting the Bellman equation {eq}`lp_bellman` when $p(dz,u)$ does not depend on $u$ shows that the optimal investment rate $x(k,u)$ **does not depend on $u$**. + +A demand shift is then a pure windfall: it tells firms nothing about future demand, so it does not change investment. + +Consequently the capital stock evolves **deterministically**, according to $k_{t+1} = k_t h(x(k_t)/k_t)$, while output is supplied inelastically and demand shocks move only prices. + +````{prf:theorem} +:label: lp_theorem2 + +Under independence, there are two possibilities for the capital stock. + +If + +```{math} +:label: lp_existence_iid +\int D(0, u) p(du) > \delta + \frac{r}{h'(\delta)} , +``` + +and $k_0 > 0$, then $k_t$ converges monotonically to the unique stationary value $k^c$ given implicitly by + +```{math} +:label: lp_kc +\int D(k^c, u) p(du) = \delta + \frac{r}{h'(\delta)} . +``` + +Otherwise capital converges monotonically to zero. +```` + +Condition {eq}`lp_kc` has a familiar interpretation. + +The left side is expected marginal revenue product of capital, which here is just the expected output price, because the marginal physical product is one. + +The right side is a **user cost of capital**: a depreciation term $\delta$ plus an interest term $r/h'(\delta)$. + +Lucas and Prescott observe that this case corresponds closely to the textbook dichotomy between short-run and long-run supply. + +In the short run, capacity is fixed and demand determines price. + +In the long run, demand fluctuations play no role at all: capacity is determined entirely by *average* demand. + +Let's verify this numerically by setting $\rho = 0$. + +```{code-cell} ipython +m_iid = create_model(ρ=0.0) +v_iid, idx_iid, k_next_iid, x_iid = solve_model(m_iid) + +# does the policy depend on u? +print("investment policy independent of u:", + np.allclose(x_iid, x_iid[:, [0]], atol=1e-10)) + +# stationary capital: where x(k) crosses δ k +def stationary_k(x_col, m): + "Capital where investment just maintains capacity." + d = x_col - m.δ * m.k + i = np.where(np.sign(d[:-1]) != np.sign(d[1:]))[0] + if len(i) == 0: + return np.nan + i = i[0] + return np.interp(0, [d[i+1], d[i]], [m.k[i+1], m.k[i]]) + +kc = stationary_k(x_iid[:, 0], m_iid) +user_cost = m_iid.δ + m_iid.r / h_prime(m_iid.δ, m_iid) +print(f"\nstationary capital k^c = {kc:.3f}") +print(f"expected price at k^c = {m_iid.a0 - m_iid.a1 * kc:.5f}") +print(f"user cost δ + r / h'(δ) = {user_cost:.5f}") +``` + +The stationary capital stock equates the expected price to the user cost of capital, as {eq}`lp_kc` requires. + +Now let's confirm that capital approaches $k^c$ monotonically, and from either direction. + +```{code-cell} ipython +def capital_path(m, idx, k0, T=60, u_index=None): + "Simulate capital, holding the demand state fixed if u_index is given." + ki = np.abs(m.k - k0).argmin() + path = np.empty(T) + j = len(m.u) // 2 if u_index is None else u_index + for t in range(T): + path[t] = m.k[ki] + ki = idx[ki, j] + return path + +fig, ax = plt.subplots(figsize=(8, 4.5)) +for k0 in (30.0, 55.0, 110.0, 150.0): + ax.plot(capital_path(m_iid, idx_iid, k0), lw=2, label=f'$k_0 = {k0:.0f}$') +ax.axhline(kc, color='k', ls='--', lw=1, label='$k^c$') +ax.set_xlabel('$t$'); ax.set_ylabel('$k_t$') +ax.legend() +plt.tight_layout() +plt.show() +``` + +Convergence is monotone and the limit does not depend on the initial capital stock. + +## Long run behavior with serially correlated demand + +The more interesting case allows demand shifts to be positively serially correlated, so that a high demand today signals high demand tomorrow. + +Now the current demand state *does* affect investment, and the capital stock is genuinely stochastic. + +To characterize the long run, Lucas and Prescott impose additional restrictions on the $\{u_t\}$ process, whose purpose is to guarantee that the distribution of $(k_t, u_t)$ settles down. + +In words, they assume that + +* from any current $u$, next period's shock lands in any non-degenerate interval with positive probability +* $u_t$ has a limiting distribution that does not depend on the initial $u_0$ and that puts positive probability on every non-degenerate interval +* $\text{Prob}\{u_{t+1} \geq x \mid u_t\}$ is strictly increasing in $u_t$, so high demand today always signals high demand tomorrow +* consumer surplus $s(k,u)$ converges uniformly as $u \to \pm\infty$ + +A Gaussian AR(1) process with $0 < \rho < 1$ satisfies these conditions, which is the example they give and the one we simulate. + +The analysis then proceeds by bounding the capital stock. + +Let $\bar v(k)$ and $\underline v(k)$ be the limits of the expected value function as $u \to \infty$ and $u \to -\infty$, and let $\bar x(k)$ and $\underline x(k)$ be the associated investment policies. + +Because investment is increasing in $u$, these bracket the policy for every $u$. + +Let $\bar k$ solve $\bar x(k) = \delta k$ and let $\underline k$ solve $\underline x(k) = \delta k$. + +These are the capital stocks that would be sustained under permanently maximal and permanently minimal demand. + +Lucas and Prescott then prove that + +* the sets $(0, \underline k)$ and $(\bar k, \infty)$, paired with any demand state, are **transient**: once the industry leaves them it does not return, and from any starting point it enters $(\underline k, \bar k)$ with probability approaching one +* the set $B = (\underline k, \bar k) \times E$ is a **single ergodic set** + +```{prf:theorem} +:label: lp_theorem3 + +If $B$ is non-empty, then for all $(k,u)$ and every initial state $(k_0, u_0)$, + +$$ +\lim_{t \to \infty} \text{Prob}\{ k_t \leq k, u_t \leq u \mid k_0, u_0 \} = P(k,u) +$$ + +exists and **does not depend on $(k_0, u_0)$**. + +The function $P$ is a probability distribution that assigns probability zero to the transient sets and positive probability to every subset of $B$ with positive area. +``` + +```{prf:theorem} +:label: lp_theorem4 + +If $B$ is non-empty, then for any initial state, with probability one + +$$ +\lim_{T \to \infty} \frac{1}{T}\sum_{t=1}^T k_t = k^* , +$$ + +where $k^*$ is the mean of $k$ under the invariant distribution $P$. +``` + +The ergodic set is non-empty -- and so this long-run behavior applies -- if and only if + +```{math} +:label: lp_existence +\lim_{u \to \infty} D(0, u) > \delta + \frac{r}{h'(\delta)} , +``` + +which says that demand must sometimes be strong enough to justify holding any capital at all. + +{prf:ref}`lp_theorem3` and {prf:ref}`lp_theorem4` are the results that make this model an object that econometricians can use. + +The first says that the model implies a well-defined stationary probability distribution for the observable time series. + +The second says that sample averages computed from a single long realization converge to the corresponding population moments of that distribution. + +Let's compute the bounds $\underline k$ and $\bar k$ for our parameterization and check that simulations behave as the theorems say. + +```{code-cell} ipython +bounds = np.array([stationary_k(x[:, j], m) for j in range(len(m.u))]) +k_lo_star, k_hi_star = bounds.min(), bounds.max() +print(f"conditional stationary capital, lowest demand state: {k_lo_star:.2f}") +print(f"conditional stationary capital, highest demand state: {k_hi_star:.2f}") +print(f"ergodic set for capital: ({k_lo_star:.2f}, {k_hi_star:.2f})") +``` + +Now we simulate the equilibrium Markov process from two very different initial capital stocks, using the *same* sequence of demand shocks. + +```{code-cell} ipython +def simulate(m, idx, k0, T=20_000, seed=0): + "Simulate the equilibrium Markov process for (k, u)." + mc = qe.MarkovChain(m.P, m.u) + u_idx = mc.simulate_indices(T, init=len(m.u)//2, random_state=seed) + ki = np.abs(m.k - k0).argmin() + k_path = np.empty(T) + for t in range(T): + k_path[t] = m.k[ki] + ki = idx[ki, u_idx[t]] + p_path = m.a0 + m.u[u_idx] - m.a1 * k_path + return k_path, p_path + +k_low, p_low = simulate(m, idx, k0=30.0, seed=1) +k_high, p_high = simulate(m, idx, k0=150.0, seed=1) + +fig, ax = plt.subplots(figsize=(9, 4.5)) +ax.plot(k_low[:250], lw=2, label='$k_0 = 30$') +ax.plot(k_high[:250], lw=2, label='$k_0 = 150$') +ax.axhline(k_lo_star, color='k', ls='--', lw=1) +ax.axhline(k_hi_star, color='k', ls='--', lw=1, label='ergodic set') +ax.set_xlabel('$t$'); ax.set_ylabel('$k_t$') +ax.legend() +plt.tight_layout() +plt.show() +``` + +Both paths are drawn into the ergodic set and then fluctuate inside it forever. + +The next figure compares the long-run distributions of capital computed from the two simulations, after discarding a burn-in sample. + +```{code-cell} ipython +burn = 2000 +fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) + +axes[0].hist(k_low[burn:], bins=60, density=True, alpha=0.5, label='from $k_0 = 30$') +axes[0].hist(k_high[burn:], bins=60, density=True, alpha=0.5, label='from $k_0 = 150$') +axes[0].set_xlabel('$k$'); axes[0].set_ylabel('density') +axes[0].set_title('invariant distribution of capital') +axes[0].legend() + +axes[1].hist(p_low[burn:], bins=60, density=True, alpha=0.5) +axes[1].set_xlabel('$p$'); axes[1].set_ylabel('density') +axes[1].set_title('invariant distribution of price') + +plt.tight_layout() +plt.show() + +print(f"mean capital from k_0 = 30: {k_low[burn:].mean():.3f}") +print(f"mean capital from k_0 = 150: {k_high[burn:].mean():.3f}") +print(f"range visited: ({k_low[burn:].min():.2f}, {k_low[burn:].max():.2f})") +``` + +The two histograms coincide, as {prf:ref}`lp_theorem3` promises, and the capital stock stays inside the ergodic set. + +Finally, here is {prf:ref}`lp_theorem4` in action: time averages from a single realization converge to the mean of the invariant distribution. + +```{code-cell} ipython +running_mean = np.cumsum(k_low) / np.arange(1, len(k_low) + 1) + +fig, ax = plt.subplots(figsize=(8, 4.5)) +ax.plot(running_mean, lw=2) +ax.axhline(k_low[burn:].mean(), color='k', ls='--', lw=1, label='$k^*$') +ax.set_xlabel('$T$'); ax.set_ylabel(r'$T^{-1}\sum_{t \leq T} k_t$') +ax.set_xscale('log') +ax.legend() +plt.tight_layout() +plt.show() +``` + +## Relation to the rational expectations lecture + +It is worth collecting the correspondences between this lecture and {doc}`rational_expectations`. + +| | {doc}`rational_expectations` | this lecture | +|---|---|---| +| uncertainty | none | Markov demand shifter $u_t$ | +| adjustment costs | quadratic, $\gamma (y'-y)^2/2$ | concave technology $k' = k h(x/k)$ | +| equilibrium object | belief $H$ with $Y' = H(Y)$ | price process $\{p_t\}$, equivalently policy $x(k,u)$ | +| equilibrium concept | $H$ is a fixed point of $\Phi$ | anticipated price distribution equals actual | +| how it is computed | planning problem, solved as an LQ problem | planning problem, solved by dynamic programming | +| what the planner maximizes | consumer plus producer surplus | discounted consumer surplus {eq}`lp_surplus` | +| equilibrium dynamics | $Y_{t+1} = \kappa_0 + \kappa_1 Y_t$ | Markov process for $(k_t, u_t)$ | +| long-run behavior | convergence to a steady state | convergence to an invariant distribution | + +The deepest common element is the strategy for computing an equilibrium. + +In both lectures, the direct approach -- guess a law of motion, compute the induced best response, and iterate -- is unreliable, because that mapping need not be a contraction. + +In both lectures, the remedy is to find a planning problem whose Euler equations coincide with the equilibrium conditions, and then to solve the planning problem by dynamic programming. + +The lecture {doc}`rational_expectations` verifies this correspondence by matching Euler equations for a particular linear-quadratic example. + +{prf:ref}`lp_theorem1` is the general statement: for this class of economies, the set of competitive equilibria and the set of solutions of the planning problem coincide, and both are singletons. + +What the baby version cannot show, because it has no uncertainty, is the payoff that Lucas and Prescott were after: an equilibrium that is a **stationary stochastic process**, with an invariant distribution and ergodic time averages. + +That is what makes it possible to confront such a model with data, and what led on to rational expectations econometrics. + +The companion lecture {doc}`optimal_growth_uncertainty` pursues exactly this theme in a one-sector growth model. + +{cite}`BrockMirman1972` prove there the counterparts of {prf:ref}`lp_theorem3` and {prf:ref}`lp_theorem4`: the distribution of capital converges to an invariant distribution that does not depend on initial conditions, and time averages along a single realization converge to population moments. + +That lecture also shows what the shadow price of capital in such a planning problem becomes in a competitive equilibrium -- namely Tobin's $q$ -- and examines a subtle question about the differentiability of the value function that the answer depends on. + +## Exercises + +```{exercise} +:label: lp_ex1 + +The user cost of capital in {eq}`lp_kc` depends on the curvature parameter $\alpha$ of the adjustment technology through $h'(\delta) = \alpha$. + +1. Explain why a **lower** $\alpha$ -- meaning stronger adjustment costs -- should reduce the long-run capital stock. +1. For the serially independent case, compute the stationary capital stock $k^c$ for $\alpha \in \{0.4, 0.6, 0.8, 1.0\}$ and verify in each case that the marginal condition {eq}`lp_kc` holds. +1. Confirm that when $\alpha = 1$ the accumulation equation is $k_{t+1} = (1-\delta)k_t + x_t$ and the user cost is the textbook $\delta + r$. +``` + +```{solution-start} lp_ex1 +:class: dropdown +``` + +A lower $\alpha$ makes $h$ more concave, so a unit of investment buys less capacity at the margin. + +Since $h'(\delta) = \alpha$, the interest component of the user cost, $r / h'(\delta) = r/\alpha$, rises as $\alpha$ falls. + +A higher user cost must be matched by a higher expected price, and since demand slopes down, that means a smaller capital stock. + +```{code-cell} ipython +print(f"{'α':>5} {'k^c':>10} {'E[price]':>12} {'user cost':>12}") +for α in (0.4, 0.6, 0.8, 1.0): + m_α = create_model(ρ=0.0, α=α, k_lo=5.0, k_hi=160.0, n_k=600) + _, _, _, x_α = solve_model(m_α) + kc_α = stationary_k(x_α[:, 0], m_α) + price = m_α.a0 - m_α.a1 * kc_α + cost = m_α.δ + m_α.r / h_prime(m_α.δ, m_α) + print(f"{α:>5.1f} {kc_α:>10.3f} {price:>12.5f} {cost:>12.5f}") +``` + +Stronger adjustment costs (lower $\alpha$) do indeed lower the long-run capital stock. + +With $\alpha = 1$ we have $h(z) = 1 - \delta + z$, so $k' = k(1 - \delta + x/k) = (1-\delta)k + x$, and $h'(\delta) = 1$, so the user cost is $\delta + r$. + +```{solution-end} +``` + +```{exercise} +:label: lp_ex2 + +{prf:ref}`lp_theorem1` says that the planner's policy **is** the competitive equilibrium. + +Verify this numerically, using the "big $K$, little $k$" logic of {doc}`rational_expectations`. + +Solve the problem of an individual price-taking firm that + +* owns capital $k_i$ and chooses $k_i'$ subject to the same accumulation technology +* takes as given the aggregate capital stock $K$, which evolves according to the planner's policy computed above +* takes as given the price $p = a_0 + u - a_1 K$, which depends on aggregate, not own, capital + +Then check that when the firm's own capital equals aggregate capital, $k_i = K$, the firm chooses exactly what the planner chooses. + +Use a coarser grid for the firm's own capital to keep the computation small. +``` + +```{solution-start} lp_ex2 +:class: dropdown +``` + +The firm's Bellman equation is + +$$ +v_i(k_i, K, u) = \max_{k_i'} \left\{ p(K,u) k_i - x(k_i, k_i') + + \beta E\left[ v_i(k_i', K', u') \mid u \right] \right\} +$$ + +where $K' $ follows the planner's law of motion. + +Note that the firm's own capital affects its revenue but not the price. + +```{code-cell} ipython +def firm_problem(m, idx_agg, n_i=80, tol=1e-8, maxit=1000, howard=20): + "Solve an individual firm's problem taking the aggregate law of motion as given." + sub = np.linspace(0, len(m.k) - 1, n_i).astype(int) # firm grid ⊂ aggregate grid + ki = m.k[sub] + n_K, n_u = len(m.k), len(m.u) + + Xi = ki[:, None] * ((ki[None, :] / ki[:, None])**(1/m.α) - (1 - m.δ)) + Ri = np.where(Xi >= 0, -Xi, -1e12) # (k_i, k_i') + price = m.a0 + m.u[None, :] - m.a1 * m.k[:, None] # (K, u) + revenue = ki[:, None, None] * price[None, :, :] # (k_i, K, u) + + v_i = np.zeros((n_i, n_K, n_u)) + u_cols = np.arange(n_u)[None, :] + for it in range(maxit): + EV = np.tensordot(v_i, m.P, axes=([2], [1])) # E[v_i(k_i', K', u') | u] + cont = EV[:, idx_agg, u_cols] # impose K' = planner's choice + obj = Ri[:, :, None, None] + m.β * cont[None, :, :, :] + pol = obj.argmax(axis=1) + v_new = revenue + np.take_along_axis(obj, pol[:, None, :, :], axis=1)[:, 0, :, :] + + for _ in range(howard): + EV = np.tensordot(v_new, m.P, axes=([2], [1])) + cont = EV[:, idx_agg, u_cols] + v_new = (revenue + Ri[np.arange(n_i)[:, None, None], pol] + + m.β * np.take_along_axis(cont, pol, axis=0)) + + if np.max(np.abs(v_new - v_i)) < tol: + v_i = v_new + break + v_i = v_new + + return ki, sub, pol + +ki, sub, pol_firm = firm_problem(m, idx) + +# compare the firm's choice with the planner's, evaluated at k_i = K +gaps = [] +for a, K_i in enumerate(sub): + for j in range(len(m.u)): + gaps.append(abs(ki[pol_firm[a, K_i, j]] - m.k[idx[K_i, j]])) +gaps = np.array(gaps) + +print(f"firm grid spacing: {np.diff(ki).mean():.3f}") +print(f"mean |firm choice - planner choice|: {gaps.mean():.3f}") +print(f"max |firm choice - planner choice|: {gaps.max():.3f}") +``` + +The discrepancies are smaller than the spacing of the firm's own capital grid. + +So the price-taking firm, responding optimally to the price process that the planner's allocation generates, chooses to do exactly what the planner does. + +That is the content of {prf:ref}`lp_theorem1`, and it is the Lucas-Prescott counterpart of the fixed point condition $H(Y) = h(Y,Y)$ in {doc}`rational_expectations`. + +```{solution-end} +``` + +```{exercise} +:label: lp_ex3 + +{prf:ref}`lp_theorem3` says that the invariant distribution does not depend on initial conditions, but it says nothing about how *wide* that distribution is. + +Investigate how serial correlation in demand affects the ergodic set. + +1. For $\rho \in \{0.0, 0.5, 0.9, 0.98\}$, compute the ergodic bounds $\underline k$ and $\bar k$. +1. Simulate each economy and compare the invariant distributions of capital. +1. Explain the pattern. Why does the $\rho = 0$ case produce a degenerate distribution for capital? +``` + +```{solution-start} lp_ex3 +:class: dropdown +``` + +Here is one solution. + +```{code-cell} ipython +fig, ax = plt.subplots(figsize=(9, 4.5)) +print(f"{'ρ':>6} {'k_lo':>9} {'k_hi':>9} {'width':>9} {'std(k)':>9}") + +for ρ in (0.0, 0.5, 0.9, 0.98): + m_ρ = create_model(ρ=ρ) + _, idx_ρ, _, x_ρ = solve_model(m_ρ) + b = np.array([stationary_k(x_ρ[:, j], m_ρ) for j in range(len(m_ρ.u))]) + k_ρ, _ = simulate(m_ρ, idx_ρ, k0=80.0, seed=3) + print(f"{ρ:>6.2f} {b.min():>9.2f} {b.max():>9.2f} " + f"{b.max()-b.min():>9.2f} {k_ρ[2000:].std():>9.3f}") + if ρ > 0: # the ρ = 0 distribution is a spike at k^c, so we omit it here + ax.hist(k_ρ[2000:], bins=50, density=True, alpha=0.45, label=f'$\\rho = {ρ}$') + +ax.set_xlabel('$k$'); ax.set_ylabel('density') +ax.legend() +plt.tight_layout() +plt.show() +``` + +The more persistent is demand, the wider the ergodic set and the more dispersed the invariant distribution of capital. + +(The figure omits $\rho = 0$, whose distribution is a spike at $k^c$ that would dwarf the others.) + +The reason is the one Lucas and Prescott emphasize. + +Investment responds to *news about future demand*, not to current demand as such. + +When $\rho = 0$, a demand shift conveys no information about the future, so investment does not respond at all, and capital converges to the single deterministic value $k^c$ of {prf:ref}`lp_theorem2`: the invariant distribution of capital is degenerate even though prices keep fluctuating. + +As $\rho$ rises, a high demand state signals a sustained period of high prices, so firms invest more, and the capital stock inherits the persistence of demand. + +```{solution-end} +``` diff --git a/lectures/optimal_growth_uncertainty.md b/lectures/optimal_growth_uncertainty.md new file mode 100644 index 000000000..5fd88ae57 --- /dev/null +++ b/lectures/optimal_growth_uncertainty.md @@ -0,0 +1,842 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: Python 3 + language: python + name: python3 +--- + +(optimal_growth_uncertainty)= +```{raw} jupyter + +``` + +# Optimal Growth Under Uncertainty and Tobin's q + +```{contents} Contents +:depth: 2 +``` + +## Overview + +This lecture studies two papers. + +The first is {cite}`BrockMirman1972`, a classic analysis of optimal one-sector growth when the production function is hit by random shocks. + +The second is {cite}`Sargent1980q`, which adds one small ingredient -- **irreversible investment** -- and uses the result to think about James Tobin's $q$ theory of investment. + +The lecture {doc}`lucas_prescott_investment` studied a competitive industry whose equilibrium is the solution of a planning problem, and whose equilibrium is a Markov process with an invariant distribution. + +Brock and Mirman ask the same questions about a one-sector growth model: + +* does the planning problem have a well-behaved solution, with continuous and monotone policy functions? +* the optimal policy makes capital a Markov process; does that process have an **invariant distribution**? +* does the distribution of capital converge to it from *any* initial condition? +* do **time averages** computed from a single realization converge to population moments? + +The last two questions are the reason both papers matter for econometrics. + +A model that answers them affirmatively delivers a stationary, ergodic stochastic process for observable time series. + +Sample moments computed from one long realization then estimate population moments, which is exactly what the rational expectations econometrics of {cite}`HansenSargent1980` requires. + +{cite}`Sargent1980q` puts the Brock-Mirman apparatus to work on a question in macroeconomics: when, and in what sense, is there an investment demand schedule relating investment to Tobin's $q$? + +Along the way, this lecture pauses over a technical issue that turns out to be central: + +* the derivative of the planner's value function with respect to capital **is** the competitive price of used capital, i.e. it is essentially Tobin's $q$ +* so the analysis stands or falls on whether that derivative exists +* irreversible investment creates a **corner**, and at a corner the standard argument for differentiability breaks down + +We'll describe the difficulty, note a gap in the published argument, and give a corrected proof. + +Let's start with some imports: + +```{code-cell} ipython +import numpy as np +import matplotlib.pyplot as plt +from collections import namedtuple +``` + +## The Brock-Mirman planning problem + +A planner chooses consumption to maximize + +```{math} +:label: bm_objective +E_0 \sum_{t=0}^\infty \beta^t u(c_t), \qquad 0 < \beta < 1, +``` + +subject to + +$$ +c_t + i_t \leq f(k_t, r_t), \qquad c_t \geq 0, \quad i_t \geq 0 , +$$ + +where $k_t$ is capital per worker, $i_t$ is investment, and $\{r_t\}$ is a sequence of independent and identically distributed shocks to production. + +Brock and Mirman assume that for each realization of the shock, $f(\cdot, r)$ is increasing, strictly concave, and satisfies the Inada conditions + +$$ +f(0, r) = 0, \qquad f'(0, r) = +\infty, \qquad f'(\infty, r) = 0 , +$$ + +and that $f$ is increasing in $r$, so that a higher shock means more output. + +The period utility function $u$ is increasing, strictly concave, and differentiable. + +In their formulation capital depreciates fully in one period, so that $k_{t+1} = i_t$ and current resources are $s_t = f(k_{t-1}, r_{t-1})$. + +The Bellman equation is + +```{math} +:label: bm_bellman +v(s) = \max_{0 \leq c \leq s} \left\{ u(c) + \beta \int v\bigl(f(s - c, r)\bigr) \, \nu(dr) \right\} , +``` + +where $\nu$ is the distribution of the shock. + +Standard arguments -- a contraction mapping, plus preservation of monotonicity and concavity -- give a unique bounded solution $v$, attained by a unique consumption policy $c = g(s)$ and investment policy $x = h(s) = s - g(s)$. + +```{prf:proposition} +:label: bm_policies + +The optimal policies $g$ and $h$ are increasing and continuous, with $g(0) = h(0) = 0$. +``` + +Brock and Mirman prove monotonicity from the first-order condition + +```{math} +:label: bm_euler +u'(g(s)) = \beta \int u'\bigl(g(f(s - g(s), r))\bigr) \, f'(s - g(s), r) \, \nu(dr) , +``` + +which is the stochastic analogue of the Euler equation. + +If $g$ did not increase with $s$, the left side of {eq}`bm_euler` would stay constant while the right side would fall, because $f$ increases and $u'$ and $f'$ decrease. + +## Long-run behavior + +The optimal policy makes capital a Markov process, + +```{math} +:label: bm_markov +k_{t+1} = h\bigl(f(k_t, r_t)\bigr) , +``` + +and Brock and Mirman devote most of their paper to the long-run behavior of this process. + +The strategy resembles the one used in {doc}`lucas_prescott_investment`. + +Consider the deterministic difference equations obtained by fixing the shock at its lowest possible value $\alpha$ and at its highest possible value $\beta_r$: + +$$ +k_{t+1} = h(f(k_t, \alpha)), \qquad k_{t+1} = h(f(k_t, \beta_r)) . +$$ + +Let $k_a$ be a positive fixed point of the first and $k_b$ a positive fixed point of the second. + +Because $f$ is increasing in the shock and $h$ is increasing, $k_a \leq k_b$. + +Brock and Mirman say that the model has a **stable configuration of fixed points** when $k_a$ and $k_b$ are positive and $k_a < k_b$. + +They then show that + +* the sets of capital stocks below $k_a$ and above $k_b$ are **transient**: the process leaves them and does not return +* the interval $[k_a, k_b]$ is where the process eventually lives + +```{prf:theorem} +:label: bm_theorem + +There is a distribution function $F$ such that the distribution $F_t$ of $k_t$ converges to $F$ uniformly, and $F$ does not depend on the initial capital stock. +``` + +{prf:ref}`bm_theorem` is the counterpart of the invariant-distribution theorem of {doc}`lucas_prescott_investment`. + +An important part of the proof is an argument about the **inverse optimal process**, which runs the Markov process backwards in time; Brock and Mirman use it to show that a stationary distribution exists and is unique without appealing to heavier machinery from probability theory. + +Convergence of distributions is not by itself enough for econometrics. + +What an econometrician needs in addition is a **mean-ergodic theorem**: a guarantee that averages computed along a single realization converge to moments of $F$, + +```{math} +:label: bm_ergodic +\frac{1}{T}\sum_{t=1}^T \phi(k_t) \to \int \phi \, dF \qquad \text{with probability one} . +``` + +For Markov processes with a unique invariant distribution and a stable configuration of fixed points, such laws of large numbers are available; {cite}`Sargent1980q` invokes the version in Doob's book to justify computing population moments of his model and comparing them with sample moments. + +This is the bridge from Brock and Mirman's theory to the rational expectations econometrics developed in the late 1970s and early 1980s. + +Without {prf:ref}`bm_theorem` and {eq}`bm_ergodic`, a likelihood function computed from a single time series would have no firm justification. + +```{note} +{cite}`BrockMirman1972` study the case in which the shocks are independent and identically distributed. + +The lecture {doc}`lucas_prescott_investment` describes a parallel set of results for a model in which the shock is serially correlated. + +In both cases the economics is the same: the state must eventually wander inside a compact set on which the process is well behaved, and the boundaries of that set are the stationary points associated with the most and least favorable shocks. +``` + +## Irreversible investment and Tobin's q + +{cite}`Sargent1980q` modifies the model in one respect: capital, once installed, cannot be eaten. + +Output can be consumed or added to the capital stock, but the conversion does not run backwards, + +```{math} +:label: q_irreversible +K_{t+1} = (1 - \delta) K_t + I_t, \qquad I_t \geq 0 , +``` + +where $0 < \delta < 1$ is the depreciation rate. + +Production is $y_t = f(K_t)\theta_t$, and the one-period utility function $u(c_t, e_t)$ is hit by a preference shock $e_t$. + +The shocks $(\theta_t, e_t)$ are independent and identically distributed over time. + +Why does this small change matter? + +In the model with reversible investment, the price of a unit of installed capital always equals the price of a unit of newly produced output, so Tobin's $q$ is identically one and firms have no investment demand schedule at all. + +The irreversibility constraint is a **friction** that lets the price of installed capital fall below the price of new capital. + +Sargent describes a competitive economy in which households own capital, rent it to firms, and trade claims to installed capital at a relative price $p_{Kt}$, which is Tobin's $q$. + +Following {cite}`Lucas_Prescott_1971` -- exactly the strategy of {doc}`lucas_prescott_investment` -- he studies the competitive equilibrium indirectly, through the planning problem that generates it. + +The planner solves + +```{math} +:label: q_bellman +v(K, \theta, e) = \max_{K' \geq (1-\delta)K} +\left\{ u\bigl(f(K)\theta + (1-\delta)K - K', e\bigr) ++ \beta E\bigl[v(K', \theta', e')\bigr] \right\} . +``` + +Let $\lambda \geq 0$ be the multiplier on the irreversibility constraint. + +The first-order condition is + +```{math} +:label: q_foc +u_c(c, e) = \beta E\bigl[v_K(K', \theta', e')\bigr] + \lambda , +\qquad \lambda \geq 0, \quad \lambda I = 0 , +``` + +and the price of installed capital is + +```{math} +:label: q_definition +q = \frac{\beta E\bigl[v_K(K', \theta', e')\bigr]}{u_c(c,e)} = 1 - \frac{\lambda}{u_c(c,e)} . +``` + +So this model implies + +* $q \leq 1$ always +* $q = 1$ exactly when investment is positive +* $q < 1$ exactly when the irreversibility constraint binds, which is when investment is zero + +Notice what {eq}`q_definition` says: **Tobin's $q$ is the derivative of the planner's value function**, normalized by marginal utility. + +Everything therefore depends on whether $v_K$ exists. + +## Differentiability of the value function + +Here is the difficulty. + +The standard tool for differentiability of a value function is the theorem of {cite}`BenvenisteScheinkman1979`. + +```{prf:theorem} Benveniste-Scheinkman +:label: bs_theorem + +Let $v$ be concave on an open set $D$ and let $K_0 \in D$. + +Suppose $W$ is concave, differentiable at $K_0$, and satisfies $W(K) \leq v(K)$ on a neighborhood of $K_0$, with $W(K_0) = v(K_0)$. + +Then $v$ is differentiable at $K_0$ and $v'(K_0) = W'(K_0)$. +``` + +The usual way to build such a $W$ is to take the optimal plan at $K_0$, **freeze** next period's capital at its optimal value $K_0'$, and let consumption absorb the change in $K$: + +```{math} +:label: q_W1 +W_1(K) = u\bigl(f(K)\theta + (1-\delta)K - K_0', e\bigr) + \beta E\bigl[v(K_0', \theta', e')\bigr] . +``` + +This is a feasible plan, so $W_1 \leq v$, with equality at $K_0$, and it is differentiable. + +But feasibility requires $K_0' \geq (1-\delta)K$, that is $K \leq K_0'/(1-\delta)$. + +When investment is positive, $K_0' > (1-\delta)K_0$ and the plan is feasible on a full neighborhood of $K_0$. + +When investment is **zero**, $K_0' = (1-\delta)K_0$ exactly, and the plan is infeasible for every $K > K_0$. + +At a corner, the standard support function is available only from the left. + +This is precisely where the published argument becomes delicate, and in our reading it has gaps. + +```{note} +Two specific difficulties. + +First, {cite}`Sargent1980q` states the derivative of the value function in the form + +$$ +v_K(K,\theta,e) = u_c(c,e)\bigl[f'(K)\theta + (1-\delta)\bigr] , +$$ + +which is obtained by substituting the first-order condition {eq}`q_foc` **with equality** into the envelope condition. + +That substitution is legitimate only where investment is positive. + +Where investment is zero, $\lambda > 0$ and the formula overstates $v_K$; the paper's own subsequent inequality records this, so the two statements are not consistent with each other. + +Second, the argument establishes differentiability of the value function by induction along the iterates $v^j = T^j v^0$, and then passes to the limit. + +That route requires knowing how the set of capital stocks at which the constraint just binds is structured. + +The published proof handles three configurations, drawn in its Figures 3-5, and appeals to an appendix that assumes the iterates are twice differentiable almost everywhere. + +The appendix in turn justifies this by differentiating the first-order condition and observing that "since the right-hand side exists almost everywhere, so does the left," which does not follow as stated, and it asserts without proof that the set of binding points has Lebesgue measure zero. +``` + +Fortunately the conclusion is true, and there is a route to it that avoids the induction on iterates entirely. + +The key observation is that at a corner a **different** feasible plan supplies the needed support function: investing nothing is feasible at *every* capital stock, so it works on both sides. + +````{prf:proposition} +:label: q_differentiability + +Assume $u$ and $f$ are increasing, strictly concave and continuously differentiable, with $u_c(0,e) = \infty$ and $f'(0) = \infty$. + +Then for each $(\theta, e)$ the value function $v(\cdot,\theta,e)$ defined by {eq}`q_bellman` is continuously differentiable on $(0,\infty)$, and + +```{math} +:label: q_envelope +v_K(K,\theta,e) = u_c(c,e) f'(K)\theta + \beta(1-\delta) E\bigl[v_K(K',\theta',e')\bigr] , +``` + +where $c$ and $K'$ are optimal at $(K,\theta,e)$. + +Where investment is positive, {eq}`q_envelope` reduces to $v_K = u_c(c,e)[f'(K)\theta + (1-\delta)]$. + +Where investment is zero, $v_K < u_c(c,e)[f'(K)\theta + (1-\delta)]$. +```` + +````{prf:proof} +**Step 0.** The operator associated with {eq}`q_bellman` maps bounded continuous concave functions into bounded continuous strictly concave functions, so $v(\cdot,\theta,e)$ is concave and continuous, and the optimal policy is single valued and continuous. + +**Step 1 (investment positive).** Suppose $I(K_0,\theta,e) > 0$, and let $K_0'$ be the optimal choice. + +Then $K_0' > (1-\delta)K_0$, so by continuity the plan that freezes next period's capital at $K_0'$ is feasible for all $K$ in a neighborhood of $K_0$. + +The function $W_1$ of {eq}`q_W1` is therefore concave, continuously differentiable, dominated by $v$, and equal to $v$ at $K_0$. + +{prf:ref}`bs_theorem` gives differentiability at $K_0$ with + +$$ +v_K(K_0,\theta,e) = u_c(c,e)\bigl[f'(K_0)\theta + (1-\delta)\bigr] . +$$ + +Because the first-order condition holds with equality here, $\beta E[v_K(K_0',\cdot)] = u_c(c,e)$, and substituting gives {eq}`q_envelope`. + +**Step 2 (investment zero, conditional on a smaller capital stock).** Suppose $I(K_0,\theta,e) = 0$, so that $K_0' = (1-\delta)K_0$ and investing nothing is optimal at $K_0$. + +Consider instead the plan that invests nothing this period and behaves optimally thereafter: + +$$ +W_2(K) = u\bigl(f(K)\theta, e\bigr) + \beta E\bigl[v\bigl((1-\delta)K, \theta', e'\bigr)\bigr] . +$$ + +Investing nothing is feasible at **every** capital stock, so $W_2 \leq v$ everywhere, with $W_2(K_0) = v(K_0)$, and $W_2$ is concave. + +Suppose that $v(\cdot,\theta',e')$ is differentiable at $(1-\delta)K_0$ for every $(\theta',e')$. + +Then $W_2$ is differentiable at $K_0$; differentiation under the expectation is legitimate because the functions $v(\cdot,\theta',e')$ are concave, hence locally Lipschitz with a common bound on a compact neighborhood. + +{prf:ref}`bs_theorem` then gives differentiability of $v$ at $K_0$, with + +$$ +v_K(K_0,\theta,e) = u_c(c,e) f'(K_0)\theta + \beta(1-\delta) E\bigl[v_K((1-\delta)K_0, \theta', e')\bigr] , +$$ + +which is {eq}`q_envelope`. + +**Step 3 (a region where the corner cannot bind).** Because $f'(0) = \infty$, there is an $\eta > 0$ such that $I(K,\theta,e) > 0$ for every $K \leq \eta$ and every $(\theta,e)$. + +The intuition is that the marginal product of capital becomes arbitrarily large as capital approaches zero, so that giving up a little consumption today buys a great deal of consumption tomorrow. + +Making this precise requires comparing the rates at which the two sides of {eq}`q_foc` diverge as $K \downarrow 0$, since with, say, logarithmic utility both sides become unbounded. + +{cite}`Sargent1980q` proves the statement in his Appendix A (his Proposition A2), under an explicit restriction on $u$ and $f$ that he imposes for this purpose. + +That result is the one ingredient of his appendix that the argument below borrows, and it is independent of the differentiability question at issue here. + +By Step 1, $v(\cdot,\theta,e)$ is therefore differentiable on $(0,\eta]$ for every $(\theta,e)$. + +**Step 4 (bootstrap).** Let + +$$ +D = \{K > 0 : v(\cdot,\theta,e) \text{ is differentiable at } K \text{ for every } (\theta,e)\} . +$$ + +We claim that if $(0, A) \subseteq D$ then $(0, A/(1-\delta)) \subseteq D$. + +Take $K < A/(1-\delta)$ and any $(\theta,e)$. + +If $I(K,\theta,e) > 0$, Step 1 applies directly. + +If $I(K,\theta,e) = 0$, then $(1-\delta)K < A$, so $v(\cdot,\theta',e')$ is differentiable at $(1-\delta)K$ for every $(\theta',e')$, and Step 2 applies. + +By Step 3 we may start the induction at $A = \eta$. + +Since $1/(1-\delta) > 1$, iterating the claim gives $(0, \eta (1-\delta)^{-n}) \subseteq D$ for every $n$, and therefore $D = (0,\infty)$. + +**Step 5 (continuity).** A concave function that is differentiable on an open interval has a continuous derivative there, so $v(\cdot,\theta,e)$ is continuously differentiable. + +Finally, where investment is zero the first-order condition holds with $\lambda > 0$, so $\beta E[v_K(K',\cdot)] < u_c$, and {eq}`q_envelope` gives $v_K < u_c[f'(K)\theta + (1-\delta)]$. +```` + +The corrected formula has a natural reading. + +Iterating {eq}`q_envelope` forward gives + +```{math} +:label: q_envelope_sum +v_K(K_t,\theta_t,e_t) = E_t \sum_{j=0}^\infty \beta^j (1-\delta)^j \, +u_c(c_{t+j}, e_{t+j}) \, f'(K_{t+j})\theta_{t+j} . +``` + +A marginal unit of capital installed today yields marginal products for as long as it survives, and each is valued at the marginal utility of consumption on the date it arrives. + +When investment is positive, the market prices that whole stream at the cost of a unit of new output, and {eq}`q_envelope_sum` collapses to $u_c[f'(K)\theta + (1-\delta)]$. + +When investment is zero, it does not, and the difference is exactly what pushes $q$ below one. + +## Computing the model + +We now solve the model of {eq}`q_bellman` with + +$$ +u(c,e) = e \ln c, \qquad f(K) = K^a . +$$ + +The shocks $\theta$ and $e$ each take two values with equal probability and are independent of each other and over time. + +One numerical detail matters a great deal. + +The irreversibility constraint says $K' \geq (1-\delta)K$, so we want $(1-\delta)K$ to be a point of the capital grid whenever $K$ is. + +Following {cite}`Sargent1980q`, we use a **geometric** grid with ratio $(1-\delta)^{1/z}$ for an integer $z$, so that $(1-\delta)K_i = K_{i-z}$ exactly. + +```{code-cell} ipython +Model = namedtuple("Model", "β δ a z K θ e W nθ ne") + +def create_model(β=0.95, δ=0.05, a=0.25, z=10, n_K=500, K_hi=20.0, + θ_vals=(0.85, 1.15), e_vals=(0.6, 1.4)): + "Geometric capital grid so that (1-δ)K is itself a grid point." + ratio = (1 - δ)**(1/z) + K = K_hi * ratio**np.arange(n_K - 1, -1, -1) + θ, e = np.array(θ_vals), np.array(e_vals) + W = np.outer(np.ones(len(θ))/len(θ), np.ones(len(e))/len(e)) + return Model(β, δ, a, z, K, θ, e, W, len(θ), len(e)) + +f = lambda m, K: K**m.a +f_prime = lambda m, K: m.a * K**(m.a - 1) +``` + +The planner chooses next period's capital from the grid, subject to positive consumption and to irreversibility. + +```{code-cell} ipython +def solve_model(m, irreversible=True, tol=1e-10, maxit=5000, howard=50): + "Value function iteration with Howard policy improvement steps." + nK = len(m.K) + C = f(m, m.K)[:, None] * m.θ[None, :] # output (K, θ) + C = C[:, None, :] + (1 - m.δ)*m.K[:, None, None] - m.K[None, :, None] + ok = C > 1e-12 # (K, K', θ) + + if irreversible: # K' ≥ (1-δ)K + irr = np.zeros((nK, nK), bool) + for i in range(nK): + irr[i, max(i - m.z, 0):] = True + ok = ok & irr[:, :, None] + + U = np.where(ok[:, :, :, None], + m.e[None, None, None, :] * np.log(np.where(ok, C, 1.0))[:, :, :, None], + -1e12) + + v = np.zeros((nK, m.nθ, m.ne)) + for it in range(maxit): + EV = np.tensordot(v, m.W, axes=([1, 2], [0, 1])) # E v(K') over (θ', e') + obj = U + m.β * EV[None, :, None, None] + idx = obj.argmax(axis=1) + v_new = np.take_along_axis(obj, idx[:, None, :, :], axis=1)[:, 0, :, :] + + for _ in range(howard): + EV = np.tensordot(v_new, m.W, axes=([1, 2], [0, 1])) + v_new = (np.take_along_axis(U, idx[:, None, :, :], axis=1)[:, 0, :, :] + + m.β * EV[idx]) + + if np.max(np.abs(v_new - v)) < tol: + v = v_new + break + v = v_new + + K_next = m.K[idx] + I = K_next - (1 - m.δ) * m.K[:, None, None] + C_pol = (f(m, m.K)[:, None, None] * m.θ[None, :, None] + + (1 - m.δ) * m.K[:, None, None] - K_next) + return v, idx, K_next, I, C_pol +``` + +To compute $q$ we need $v_K$, and {prf:ref}`q_differentiability` tells us how to get it: equation {eq}`q_envelope` is a linear fixed point problem in $v_K$ given the optimal policy, and it is a contraction with modulus $\beta(1-\delta)$. + +```{code-cell} ipython +def marginal_value(m, idx, C_pol, tol=1e-13, maxit=50_000): + "Solve the envelope equation v_K = u_c f'(K)θ + β(1-δ) E v_K(K')." + vK = np.zeros_like(C_pol) + direct = (m.e[None, None, :]/C_pol) * f_prime(m, m.K)[:, None, None] * m.θ[None, :, None] + for _ in range(maxit): + EvK = np.tensordot(vK, m.W, axes=([1, 2], [0, 1])) + new = direct + m.β * (1 - m.δ) * EvK[idx] + if np.max(np.abs(new - vK)) < tol: + return new + vK = new + return vK + +m = create_model() +v, idx, K_next, I, C_pol = solve_model(m) +vK = marginal_value(m, idx, C_pol) + +u_c = m.e[None, None, :] / C_pol +E_vK = np.tensordot(vK, m.W, axes=([1, 2], [0, 1])) +q = m.β * E_vK[idx] / u_c +corner = I <= 1e-12 + +print(f"share of states with zero investment: {corner.mean():.3f}") +``` + +Let's check {prf:ref}`q_differentiability` numerically, by comparing the envelope formula with a finite-difference derivative of the computed value function. + +```{code-cell} ipython +vK_fd = np.gradient(v, m.K, axis=0) +naive = u_c * (f_prime(m, m.K)[:, None, None] * m.θ[None, :, None] + (1 - m.δ)) +mid = ((m.K > 2.2) & (m.K < 7.0))[:, None, None] & np.ones_like(corner) + +rel = lambda x: np.median((np.abs(x - vK_fd)/np.abs(vK_fd))[mid & sel]) + +sel = ~corner +print(f"investment positive: envelope error {rel(vK):.4f}, " + f"naive formula error {rel(naive):.4f}") +sel = corner +print(f"investment zero: envelope error {rel(vK):.4f}, " + f"naive formula error {rel(naive):.4f}") +``` + +Where investment is positive, both formulas agree with the numerical derivative. + +Where investment is zero, the envelope formula {eq}`q_envelope` remains accurate while the formula $u_c[f'(K)\theta + (1-\delta)]$ is off by tens of percent. + +This is the practical content of the correction: at corners the two expressions are genuinely different objects, and it is the envelope formula that gives the price of installed capital. + +We can also verify Step 3 of the proof, which asserts a region of small capital stocks where the corner never binds. + +```{code-cell} ipython +m_wide = create_model(n_K=700) # a grid reaching lower capital stocks +_, idx_w, _, I_w, _ = solve_model(m_wide) +frac_corner = (I_w <= 1e-12).reshape(len(m_wide.K), -1).mean(axis=1) +first = np.argmax(frac_corner > 0) +print(f"grid runs from K = {m_wide.K.min():.3f}") +print(f"investment is positive for every shock when K < {m_wide.K[first]:.3f}") +``` + +### Policies and the price of installed capital + +```{code-cell} ipython +fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) +show = (m.K > 1.5) & (m.K < 9) + +for j in range(m.nθ): + for l in range(m.ne): + lab = rf'$\theta={m.θ[j]}, e={m.e[l]}$' + axes[0].plot(m.K[show], K_next[show, j, l], label=lab) + axes[1].plot(m.K[show], q[show, j, l], label=lab) + +axes[0].plot(m.K[show], m.K[show], 'k--', lw=1, label='45 degree line') +axes[0].plot(m.K[show], (1-m.δ)*m.K[show], 'k:', lw=1, label=r'$(1-\delta)K$') +axes[0].set_xlabel('$K$'); axes[0].set_ylabel("$K'$"); axes[0].set_title('law of motion') +axes[1].axhline(1.0, color='k', ls='--', lw=1) +axes[1].set_xlabel('$K$'); axes[1].set_ylabel('$q$'); axes[1].set_title("Tobin's $q$") +axes[0].legend(fontsize=8); axes[1].legend(fontsize=8) +plt.tight_layout() +plt.show() +``` + +The law of motion runs along the lower dotted line $(1-\delta)K$ whenever the irreversibility constraint binds. + +Exactly on that region, $q$ falls below one. + +### The invariant distribution + +As in {doc}`lucas_prescott_investment`, we check that the long-run distribution does not depend on where the economy starts. + +```{code-cell} ipython +def simulate(m, idx, K0, T=50_000, seed=0, burn=500): + "Simulate the equilibrium Markov process." + rng = np.random.default_rng(seed) + j_idx = rng.integers(0, m.nθ, T) + l_idx = rng.integers(0, m.ne, T) + k = np.abs(m.K - K0).argmin() + path = np.empty(T, int) + for t in range(T): + path[t] = k + k = idx[k, j_idx[t], l_idx[t]] + return path[burn:], j_idx[burn:], l_idx[burn:] + +p_lo, j_lo, l_lo = simulate(m, idx, K0=2.0, seed=7) +p_hi, j_hi, l_hi = simulate(m, idx, K0=15.0, seed=7) + +print(f"mean capital starting from K0 = 2: {m.K[p_lo].mean():.4f}") +print(f"mean capital starting from K0 = 15: {m.K[p_hi].mean():.4f}") +print(f"capital visited: [{m.K[p_lo].min():.2f}, {m.K[p_lo].max():.2f}]") +``` + +```{code-cell} ipython +fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) + +axes[0].hist(m.K[p_lo], bins=60, density=True, alpha=0.5, label='from $K_0 = 2$') +axes[0].hist(m.K[p_hi], bins=60, density=True, alpha=0.5, label='from $K_0 = 15$') +axes[0].set_xlabel('$K$'); axes[0].set_ylabel('density') +axes[0].set_title('invariant distribution of capital') +axes[0].legend() + +q_sim, I_sim = q[p_lo, j_lo, l_lo], I[p_lo, j_lo, l_lo] +axes[1].scatter(q_sim, I_sim, s=4, alpha=0.2) +axes[1].set_xlabel('$q$'); axes[1].set_ylabel('$I$') +axes[1].set_title("investment and Tobin's $q$") + +plt.tight_layout() +plt.show() +``` + +The two histograms coincide, as {prf:ref}`bm_theorem` leads us to expect. + +The scatter plot on the right reproduces the central picture of {cite}`Sargent1980q`. + +Investment is positive only when $q$ equals one; when $q$ is below one, investment is exactly zero. + +```{code-cell} ipython +print(f"q when I > 0: [{q_sim[I_sim > 1e-12].min():.4f}, {q_sim[I_sim > 1e-12].max():.4f}]") +print(f"q when I = 0: [{q_sim[I_sim <= 1e-12].min():.4f}, {q_sim[I_sim <= 1e-12].max():.4f}]") +``` + +The small departures from exactly one reflect the discreteness of the capital grid. + +Sargent's point is now visible. + +There *is* a positive relationship between investment and $q$ in these data, and an econometrician could certainly run that regression. + +But the relationship is not an investment demand schedule. + +It is a **mongrel relation** that mixes together preferences, technology and the distribution of the shocks, and it will shift whenever any of them changes. + +{ref}`ogu_ex3` asks you to demonstrate this. + +## Exercises + +```{exercise} +:label: ogu_ex1 + +{cite}`BrockMirman1972` is famous partly for a special case that can be solved by hand. + +Let $u(c) = \ln c$, let $f(k, r) = r k^\alpha$ with $0 < \alpha < 1$, and let capital depreciate fully each period, so that $k_{t+1} = f(k_t, r_t) - c_t$. + +1. Verify that the optimal policy is $k_{t+1} = \alpha\beta\, r_t k_t^\alpha$. +1. Show that $\ln k_t$ follows an AR(1) process, and derive the mean and variance of its invariant distribution when $\ln r_t \sim N(\mu, \sigma^2)$. +1. Simulate the model and check both the invariant distribution and the mean-ergodic property {eq}`bm_ergodic`. +``` + +```{solution-start} ogu_ex1 +:class: dropdown +``` + +Guess that a constant fraction of resources is saved, $k_{t+1} = s\, r_t k_t^\alpha$. + +The Euler equation for this problem is + +$$ +\frac{1}{c_t} = \beta E_t \left[ \frac{1}{c_{t+1}} \alpha r_{t+1} k_{t+1}^{\alpha - 1} \right] . +$$ + +With $c_t = (1-s) r_t k_t^\alpha$ and $k_{t+1} = s r_t k_t^\alpha$, the right side becomes + +$$ +\beta E_t \left[\frac{\alpha r_{t+1}k_{t+1}^{\alpha-1}}{(1-s) r_{t+1} k_{t+1}^{\alpha}}\right] += \frac{\beta \alpha}{(1-s) k_{t+1}} += \frac{\beta\alpha}{(1-s) s r_t k_t^{\alpha}} , +$$ + +while the left side is $1/[(1-s) r_t k_t^\alpha]$. + +Equating the two gives $s = \alpha\beta$. + +Taking logs of the policy, + +$$ +\ln k_{t+1} = \ln(\alpha\beta) + \alpha \ln k_t + \ln r_t , +$$ + +an AR(1) with coefficient $\alpha$. + +With $\ln r_t \sim N(\mu, \sigma^2)$, the invariant distribution of $\ln k$ is normal with + +$$ +\text{mean} = \frac{\ln(\alpha\beta) + \mu}{1 - \alpha}, +\qquad +\text{variance} = \frac{\sigma^2}{1 - \alpha^2} . +$$ + +```{code-cell} ipython +α, β_d, μ, σ = 0.4, 0.95, 0.0, 0.1 +T = 200_000 +rng = np.random.default_rng(0) +ln_r = μ + σ * rng.normal(size=T) + +k = np.empty(T + 1) +k[0] = 1.0 +for t in range(T): + k[t+1] = α * β_d * np.exp(ln_r[t]) * k[t]**α + +ln_k = np.log(k[1000:]) +print(f"mean of ln k: simulated {ln_k.mean():.4f}, " + f"theory {(np.log(α*β_d) + μ)/(1-α):.4f}") +print(f"var of ln k: simulated {ln_k.var():.4f}, " + f"theory {σ**2/(1-α**2):.4f}") +``` + +The time averages match the population moments of the invariant distribution, which is the mean-ergodic property {eq}`bm_ergodic` in action. + +```{solution-end} +``` + +```{exercise} +:label: ogu_ex2 + +This exercise examines the difference between the two candidate formulas for $v_K$ discussed above. + +For the baseline model, compute at every state + +* the envelope formula {eq}`q_envelope` +* the formula $u_c(c,e)[f'(K)\theta + (1-\delta)]$ +* a finite-difference derivative of the value function + +Then + +1. Confirm that all three agree where investment is positive. +1. Confirm that the second disagrees with the other two where investment is zero, and report by how much. +1. Explain, using {eq}`q_foc`, why the second formula is an upper bound for $v_K$. +1. What would go wrong with the computed $q$ if you used the second formula everywhere? +``` + +```{solution-start} ogu_ex2 +:class: dropdown +``` + +The second formula is what you get by substituting the first-order condition {eq}`q_foc` into the envelope condition **as though** $\lambda = 0$. + +Since $\lambda \geq 0$, dropping it can only raise the expression, so + +$$ +v_K = u_c f'(K)\theta + (1-\delta)\bigl[u_c - \lambda\bigr] +\leq u_c\bigl[f'(K)\theta + (1-\delta)\bigr] , +$$ + +with equality exactly when $\lambda = 0$, that is, when investment is positive. + +```{code-cell} ipython +gap = np.abs(naive - vK_fd)/np.abs(vK_fd) +env = np.abs(vK - vK_fd)/np.abs(vK_fd) + +for name, sel in (("investment positive", ~corner), ("investment zero", corner)): + s = mid & sel + print(f"{name:20}: envelope {np.median(env[s]):.4f} " + f"naive {np.median(gap[s]):.4f} max naive {gap[s].max():.4f}") + +q_naive = m.β * np.tensordot(naive, m.W, axes=([1, 2], [0, 1]))[idx] / u_c +print(f"\nusing the naive formula, max q = {q_naive.max():.3f} " + f"(theory says q ≤ 1)") +print(f"fraction of states with q > 1: {(q_naive > 1 + 1e-8).mean():.3f}") +``` + +Using the second formula everywhere produces values of $q$ that exceed one, which the theory rules out. + +The reason is economic, not numerical: at a corner the market does not price installed capital at the cost of new capital, precisely because the household would like to sell capital and cannot. + +```{solution-end} +``` + +```{exercise} +:label: ogu_ex3 + +Show that the regression of investment on $q$ is not structural. + +Solve and simulate the model for several economies that differ in + +* the amount of preference risk (the spread of $e$) +* the amount of technology risk (the spread of $\theta$) +* the depreciation rate $\delta$ + +For each economy, report the slope and $R^2$ of a regression of $I$ on $q$, together with the frequency with which investment is zero. + +Comment on what this implies for an econometrician who estimates an "investment demand schedule" relating investment to $q$. +``` + +```{solution-start} ogu_ex3 +:class: dropdown +``` + +Here is one solution. + +```{code-cell} ipython +def q_regression(**kwargs): + "Solve an economy, simulate it, and regress I on q." + mm = create_model(**kwargs) + vv, ii, _, II, CC = solve_model(mm) + vvK = marginal_value(mm, ii, CC) + uu_c = mm.e[None, None, :]/CC + qq = mm.β * np.tensordot(vvK, mm.W, axes=([1, 2], [0, 1]))[ii] / uu_c + p, j, l = simulate(mm, ii, K0=3.5) + q_s, I_s = qq[p, j, l], II[p, j, l] + slope = np.polyfit(q_s, I_s, 1)[0] + corr = np.corrcoef(I_s, q_s)[0, 1] + return slope, corr**2, (I_s <= 1e-12).mean() + +cases = [("baseline", {}), + ("less preference risk", dict(e_vals=(0.8, 1.2))), + ("more preference risk", dict(e_vals=(0.4, 1.6))), + ("less technology risk", dict(θ_vals=(0.95, 1.05))), + ("more technology risk", dict(θ_vals=(0.7, 1.3))), + ("faster depreciation", dict(δ=0.10))] + +print(f"{'economy':24}{'slope':>9}{'R^2':>8}{'zero I':>9}") +for name, kw in cases: + slope, r2, zero = q_regression(**kw) + print(f"{name:24}{slope:9.3f}{r2:8.3f}{zero:9.3f}") +``` + +The slope and the fit both move around substantially across economies that share the same preferences and technology parameters but differ in the distribution of shocks or in the depreciation rate. + +Nothing in the regression is invariant, so it cannot be used to predict what investment would do under a policy that changed any of these features. + +This is the message of {cite}`Sargent1980q`, and it is a concrete instance of the Lucas critique: the very friction that makes $q$ interesting -- the occasionally binding irreversibility constraint -- also makes the relationship between $I$ and $q$ a reduced-form artifact rather than a decision rule. + +```{solution-end} +``` diff --git a/lectures/rational_expectations.md b/lectures/rational_expectations.md index 844d98286..1665fe366 100644 --- a/lectures/rational_expectations.md +++ b/lectures/rational_expectations.md @@ -46,6 +46,14 @@ due to Lucas and Prescott {cite}`Lucas_Prescott_1971`. That 1971 paper is one of a small number of research articles that ignited a *rational expectations revolution*. +You can regard the present lecture as a "baby" version of that paper, and as an introduction to ideas that two sequels study in more general settings. + +{doc}`lucas_prescott_investment` describes the model that Lucas and Prescott actually built, in which demand is shifted by a Markov process, so that the equilibrium is a stochastic process rather than a deterministic path. + +{doc}`optimal_growth_uncertainty` studies a one-sector optimal growth model with random shocks to production, and uses it to think about Tobin's $q$. + +Both sequels ask a question that cannot arise here, because this lecture has no uncertainty: does the equilibrium have an invariant probability distribution to which it converges from any initial condition? + We follow Lucas and Prescott by employing a setting that is readily "Bellmanized" (i.e., susceptible to being formulated as a dynamic programming problems. Because we use linear quadratic setups for demand and costs, we can deploy the LQ programming techniques described in {doc}`this lecture `. @@ -207,6 +215,8 @@ References for this lecture include * {cite}`Sargent1987`, chapter XIV * {cite}`Ljungqvist2012`, chapter 7 +The two sequels to this lecture are {doc}`lucas_prescott_investment` and {doc}`optimal_growth_uncertainty`. + ## Rational expectations equilibrium ```{index} single: Rational Expectations Equilibrium; Definition @@ -578,6 +588,30 @@ y_{t+1} = h_0 + h_1 y_t + h_2 Y_t Hence a rational expectations equilibrium will be defined by the parameters $(\kappa_0, \kappa_1, h_0, h_1, h_2)$ in {eq}`ree_hlom2`--{eq}`ree_ex5`. +## Concluding remarks + +Three ideas from this lecture recur throughout macroeconomics. + +The first is the equilibrium concept itself: a **perceived** law of motion for a market-wide object must coincide with the **actual** law of motion that the resulting decisions generate. + +The second is the "Big $Y$, little $y$" device that lets a price-taking firm be representative. + +The third is the computational strategy: because the mapping from perceived to actual laws of motion is not a contraction, we found an equilibrium by solving a **planning problem** instead, and then read off equilibrium prices as shadow prices. + +The two sequels take these ideas into settings with uncertainty. + +{doc}`lucas_prescott_investment` returns to {cite}`Lucas_Prescott_1971` itself. + +There, demand is shifted by a Markov process, firms face a nonlinear technology for converting investment into capacity, and the equilibrium is proved to exist, to be unique, and to solve a planning problem that maximizes discounted consumer surplus. + +Because the equilibrium is a Markov process, one can ask whether it converges to an invariant distribution, and whether time averages computed from a single realization converge to population moments. + +Affirmative answers are what make such models usable in econometrics. + +{doc}`optimal_growth_uncertainty` pursues the same questions in the one-sector optimal growth model of {cite}`BrockMirman1972`, and then follows {cite}`Sargent1980q` in adding irreversible investment. + +That small change makes the shadow price of installed capital -- Tobin's $q$ -- diverge from the price of new capital, and it turns the planner's value function into the object on which a theory of investment rests. + ## Exercises ```{exercise} diff --git a/lectures/two_auctions.md b/lectures/two_auctions.md index dd8c40215..3ba9688e8 100644 --- a/lectures/two_auctions.md +++ b/lectures/two_auctions.md @@ -13,7 +13,7 @@ kernelspec: # First-Price and Second-Price Auctions -This lecture is designed to set the stage for a subsequent lecture about [Multiple Good Allocation Mechanisms](https://python.quantecon.org/house_auction.html) +This lecture is designed to set the stage for a subsequent lecture about {doc}`house_auction` In that lecture, a planner or auctioneer simultaneously allocates several goods to set of people. @@ -63,7 +63,7 @@ Much of our Python code below is based on his. **Detailed Setting:** -There are $n>2$ prospective buyers named $i = 1, 2, \ldots, n$. +There are $n \geq 2$ prospective buyers named $i = 1, 2, \ldots, n$. Buyer $i$ attaches value $v_i$ to the good being sold. @@ -93,7 +93,15 @@ Bidder optimally chooses to bid less than $v_i$. ### Characterization of FPSB auction -A FPSB auction has a unique symmetric Bayesian Nash Equilibrium. +We assume throughout that + +* valuations are **private** and **independent** across bidders +* bidders are **symmetric**: their valuations are drawn from a common distribution $F$ that is continuous and strictly increasing on its support +* bidders are **risk neutral** + +Under these assumptions a FPSB auction has a unique Bayesian Nash equilibrium in symmetric, strictly increasing bidding strategies. + +Because the equilibrium bidding strategy is strictly increasing, the bidder with the highest valuation submits the highest bid, and so wins. The optimal bid of buyer $i$ is @@ -109,7 +117,9 @@ $$ (eq:optbid2) -A proof for this assertion is available at the [Wikipedia page](https://en.wikipedia.org/wiki/Vickrey_auction) about Vickrey auctions +For a derivation, see the [Wikipedia page](https://en.wikipedia.org/wiki/First-price_sealed-bid_auction) about first-price sealed-bid auctions, or {cite}`Krishna2009`, chapter 2. + +We'll verify this formula by simulation below, and {ref}`ta_ex2` asks you to derive an equivalent expression that is easy to evaluate for any distribution $F$. +++ @@ -123,7 +133,11 @@ A proof for this assertion is available at the [Wikipedia page](https://en.wiki In a SPSB auction bidders optimally choose to bid their values. -Formally, a dominant strategy profile in a SPSB auction with a single, indivisible item has each bidder bidding its value. +Formally, in a SPSB auction with a single, indivisible item, bidding one's own value is a **weakly dominant** strategy. + +It is *weakly* dominant because a bidder who bids something other than her value never does better, and sometimes does worse, whatever the other bidders do. + +Notice how much stronger this is than the FPSB result: it requires no assumption at all about the distribution of other bidders' valuations, nor about how they bid. A proof is provided at [the Wikipedia page](https://en.wikipedia.org/wiki/Vickrey_auction) about Vickrey auctions @@ -142,7 +156,17 @@ We'll simulate outcomes and, by using a law of large numbers, verify that the We can use our simulation to illustrate a **Revenue Equivalence Theorem** that asserts that on average first-price and second-price sealed bid auctions provide a seller the same revenue. -To read about the revenue equivalence theorem, see [this Wikepedia page](https://en.wikipedia.org/wiki/Revenue_equivalence) +The theorem requires hypotheses that both of our auctions satisfy: + +* valuations are independent and private, and bidders are symmetric and risk neutral +* the two mechanisms award the good to the bidder with the highest valuation +* a bidder with the lowest possible valuation expects zero surplus + +Under these hypotheses any two such mechanisms yield the same expected payment for each bidder, and therefore the same expected revenue for the seller. + +{ref}`ta_ex4` shows what happens when one of these hypotheses -- risk neutrality -- fails. + +To read about the revenue equivalence theorem, see [this Wikipedia page](https://en.wikipedia.org/wiki/Revenue_equivalence) +++ @@ -197,7 +221,7 @@ import scipy.stats as stats import scipy.interpolate as interp # for plots -plt.rcParams.update({"text.usetex": True, 'font.size': 14}) +plt.rcParams.update({'font.size': 14}) colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] # ensure the notebook generates the same randomness @@ -231,13 +255,8 @@ idx = np.argsort(v, axis=0) v = np.take_along_axis(v, idx, axis=0) b = np.take_along_axis(b, idx, axis=0) -# the id for the bidders is created. -ii = np.repeat(np.arange(1, N+1)[:, None], R, axis=1) -# the id is sorted according to bid price as well. -ii = np.take_along_axis(ii, idx, axis=0) - -# In FPSB and SPSB, winners are those with highest values. -winning_player = ii[-1, :] +# In FPSB and SPSB the winner is the bidder with the highest valuation, +# which after sorting is the last row. # highest bid winner_pays_fpsb = b[-1, :] @@ -389,21 +408,27 @@ def evaluate_largest(v_hat, array, order=1): A method to estimate the largest (or certain-order largest) value of the other biders, conditional on player 1 wins the auction. + We estimate E[y | y < v_hat], where y is the highest valuation among the + other bidders. We do this by taking bidder 1 as the reference bidder + (valuations are i.i.d., so the choice does not matter), discarding her row, + and averaging the highest remaining valuation over those auctions in which + every other bidder's valuation falls below v_hat. + Parameters: ---------- - v_hat : float, the value of player 1. The biggest value in the auction that player 1 wins. + v_hat : float, the valuation of the reference bidder. array: 2 dimensional array of bidders' values in shape of (N,R), where N: number of players, R: number of auctions - order: int. The order of largest number among bidders who lose. - e.g. the order for largest number beside winner is 1. - the order for second-largest number beside winner is 2. + order: int. Which order statistic of the losing bidders to average. + order=1 gives the highest losing valuation, + order=2 the second highest, and so on. """ N, R = array.shape - # drop the first row because we assume first row is the winner's bid + # discard the reference bidder's row; condition on the rest losing array_residual = array[1:, :].copy() winning_auctions_mask = (array_residual < v_hat).all(axis=0) @@ -511,11 +536,6 @@ idx = np.argsort(v, axis=0) v = np.take_along_axis(v, idx, axis=0) b = np.take_along_axis(b, idx, axis=0) -ii = np.repeat(np.arange(1, N + 1)[:, None], R, axis=1) -ii = np.take_along_axis(ii, idx, axis=0) - -winning_player = ii[-1, :] - # highest bid winner_pays_fpsb = b[-1, :] # 2nd-highest valuation @@ -592,7 +612,7 @@ class bid_price_solution: return np.mean(order_largest_bids) - def compute_optimal_bid_FPSB(self): + def compute_optimal_bid_FPSB(self, plot=True): # we compute the quantile of v as our grid pct_quantile = np.linspace(0, 100, 101)[1:-1] v_grid = np.percentile(self.value_mat.flatten(), q=pct_quantile) @@ -607,6 +627,9 @@ class bid_price_solution: self.b_star_num = interp.interp1d(v_grid, EV, fill_value="extrapolate") + if not plot: + return None + pct_quantile_fine = np.linspace(0, 100, 1001)[1:-1] v_grid_fine = np.percentile(self.value_mat.flatten(), q=pct_quantile_fine) @@ -625,6 +648,9 @@ class bid_price_solution: return None def plot_winner_payment_distribution(self): + if not hasattr(self, 'b_star_num'): # bids have not been computed yet + self.compute_optimal_bid_FPSB(plot=False) + self.b = self.b_star_num(self.value_mat) idx = np.argsort(self.value_mat, axis=0) @@ -632,12 +658,6 @@ class bid_price_solution: self.v = np.take_along_axis(self.value_mat, idx, axis=0) self.b = np.take_along_axis(self.b, idx, axis=0) - N, R = self.value_mat.shape - self.ii = np.repeat(np.arange(1, N + 1)[:, None], R, axis=1) - self.ii = np.take_along_axis(self.ii, idx, axis=0) - - winning_player = self.ii[-1, :] - # highest bid winner_pays_fpsb = self.b[-1, :] # 2nd-highest valuation @@ -681,13 +701,246 @@ chi_squ_case.compute_optimal_bid_FPSB() chi_squ_case.plot_winner_payment_distribution() ``` -## References +## Exercises -+++ +```{exercise} +:label: ta_ex1 + +Verify the Revenue Equivalence Theorem by simulation. + +For $n = 2, 3, 5, 10$ bidders with valuations drawn independently from $U(0,1)$, simulate many auctions and compute + +1. the average payment of the winner in a FPSB auction, in which each bidder bids $\frac{n-1}{n} v_i$ +1. the average payment of the winner in a SPSB auction, in which each bidder bids $v_i$ + +Compare both with the theoretical expected revenue $\frac{n-1}{n+1}$, and comment on how the seller's revenue changes with the number of bidders. +``` + +```{solution-start} ta_ex1 +:class: dropdown +``` + +```{code-cell} ipython3 +R_ex = 200_000 +rng_ex = np.random.default_rng(1234) + +print(f"{'n':>4}{'FPSB':>12}{'SPSB':>12}{'(n-1)/(n+1)':>14}") +for n in (2, 3, 5, 10): + v_ex = np.sort(rng_ex.uniform(0, 1, (n, R_ex)), axis=0) + fpsb = (n - 1)/n * v_ex[-1, :] # winner's own bid + spsb = v_ex[-2, :] # second highest valuation + print(f"{n:>4}{fpsb.mean():>12.4f}{spsb.mean():>12.4f}{(n-1)/(n+1):>14.4f}") +``` + +The two auctions raise the same expected revenue, and both converge to the highest possible valuation as $n$ grows. + +With more bidders, competition pushes the winning payment toward the top of the support of valuations. + +Notice that the two auctions raise the same revenue on average even though the *distributions* of the winner's payment differ: in a FPSB auction the payment is a deterministic function of the winner's valuation, while in a SPSB auction it is the second-highest valuation, which is random given the winner's valuation. + +```{solution-end} +``` + +```{exercise} +:label: ta_ex2 + +Equation {eq}`eq:optbid1` says that an optimal bid in a FPSB auction is $\mathbf{E}[y_i \mid y_i < v_i]$. + +1. Show that this can be written + + $$ + b(v) = v - \frac{\int_0^{v} F(x)^{n-1} dx}{F(v)^{n-1}} + $$ + + where $F$ is the distribution function of a valuation. + +1. Verify that this reduces to $\frac{n-1}{n}v$ when $F$ is uniform on $[0,1]$. + +1. Evaluate the formula for a $\chi^2(2)$ distribution of valuations and compare it with the simulation-based bid function computed in the lecture. +``` + +```{solution-start} ta_ex2 +:class: dropdown +``` + +The distribution function of $y_i = \max_{j \neq i} v_j$ is $\tilde F_{n-1}(y) = F(y)^{n-1}$. + +Hence + +$$ +\mathbf{E}[y \mid y < v] = \frac{1}{F(v)^{n-1}} \int_0^v y \, d\left[F(y)^{n-1}\right] . +$$ + +Integrating by parts, + +$$ +\int_0^v y \, d\left[F(y)^{n-1}\right] = v F(v)^{n-1} - \int_0^v F(y)^{n-1} dy , +$$ + +which gives the formula. + +For $F(x) = x$ on $[0,1]$ we get $b(v) = v - \frac{v^n/n}{v^{n-1}} = \frac{n-1}{n} v$. + +The formula has a nice reading: a bidder shades her bid below her valuation by an amount that shrinks as the number of competitors grows. + +```{code-cell} ipython3 +from scipy.integrate import quad + +def b_closed_form(v, F, n): + "Optimal FPSB bid for a bidder with valuation v when rivals' values ~ F." + shading = quad(lambda x: F(x)**(n - 1), 0, v)[0] / F(v)**(n - 1) + return v - shading + +# check against the analytical solution for the uniform case +print("uniform check") +for v0 in (0.3, 0.6, 0.9): + print(f" v = {v0}: closed form {b_closed_form(v0, lambda x: x, N):.4f}, " + f"analytical {b_star(v0, N):.4f}") +``` + +```{code-cell} ipython3 +# now the chi-squared case studied in the lecture +F_chi2 = stats.chi2(df=2).cdf +v_test = np.percentile(v.flatten(), [10, 30, 50, 70, 90]) + +print(f"{'v':>8}{'closed form':>14}{'simulated':>12}") +for v0 in v_test: + print(f"{v0:>8.3f}{b_closed_form(v0, F_chi2, N):>14.4f}" + f"{float(b_star_num(v0)):>12.4f}") +``` + +The closed form and the simulation agree closely, which is a useful check on both. + +```{solution-end} +``` + +```{exercise} +:label: ta_ex3 + +This exercise asks you to see *why* truthful bidding is a weakly dominant strategy in a SPSB auction but not in a FPSB auction. + +Fix $n = 5$ and consider a bidder whose valuation is $v = 0.75$, facing rivals whose valuations are $U(0,1)$. + +1. In a SPSB auction the rivals bid truthfully. Compute this bidder's expected surplus as a function of her own bid $b$ and plot it. +1. In a FPSB auction the rivals bid $\frac{n-1}{n}v_j$. Compute and plot her expected surplus as a function of $b$. +1. Where does each curve peak? What surplus does she earn in the FPSB auction if she bids her valuation? +``` + +```{solution-start} ta_ex3 +:class: dropdown +``` + +In the SPSB auction she wins when $y < b$ and then pays $y$, so her expected surplus is + +$$ +\int_0^b (v - y) \, (n-1) y^{n-2} dy . +$$ + +Differentiating with respect to $b$ gives $(v-b)(n-1)b^{n-2}$, which is positive for $b < v$ and negative for $b > v$, so $b = v$ is optimal. + +In the FPSB auction she wins when every rival's bid falls below $b$, which happens with probability $\left(\frac{nb}{n-1}\right)^{n-1}$, and then she pays $b$. + +```{code-cell} ipython3 +n_ex, v_own = 5, 0.75 +bids = np.linspace(0, 1, 401) + +spsb_surplus = [quad(lambda y: (v_own - y)*(n_ex - 1)*y**(n_ex - 2), + 0, min(bb, 1))[0] for bb in bids] +fpsb_surplus = [(v_own - bb)*min(1, n_ex*bb/(n_ex - 1))**(n_ex - 1) + for bb in bids] + +fig, ax = plt.subplots(figsize=(6, 4)) +ax.plot(bids, spsb_surplus, label='SPSB') +ax.plot(bids, fpsb_surplus, label='FPSB') +ax.axvline(v_own, ls='--', c='k', lw=1, label='own valuation') +ax.axvline((n_ex - 1)/n_ex*v_own, ls=':', c='r', lw=1, label='FPSB optimal bid') +ax.set_xlabel('own bid $b$') +ax.set_ylabel('expected surplus') +ax.legend() +plt.show() + +print(f"SPSB surplus is maximized at b = {bids[int(np.argmax(spsb_surplus))]:.3f}") +print(f"FPSB surplus is maximized at b = {bids[int(np.argmax(fpsb_surplus))]:.3f}" + f" (theory: {(n_ex - 1)/n_ex*v_own:.3f})") +print(f"FPSB surplus from bidding one's valuation: {(v_own - v_own):.3f}") +``` + +The SPSB curve peaks exactly at the bidder's valuation. + +The FPSB curve peaks strictly below it, and bidding one's valuation in a FPSB auction earns a surplus of exactly zero: the bidder wins more often, but pays her full valuation whenever she does. + +```{solution-end} +``` + +```{exercise} +:label: ta_ex4 + +The Revenue Equivalence Theorem requires bidders to be **risk neutral**. + +Suppose instead that each bidder has utility $u(x) = x^\rho$ with $0 < \rho \leq 1$, so that $\rho < 1$ means risk aversion, and that valuations are $U(0,1)$. + +One can show that the symmetric equilibrium bid in a FPSB auction becomes + +$$ +b(v) = \frac{n-1}{n-1+\rho} v . +$$ + +1. Verify this numerically: for $n = 5$ and a bidder with $v = 0.8$, compute expected utility as a function of her own bid when rivals use this rule, and check where it is maximized. +1. Compute the seller's expected revenue in FPSB and SPSB for $\rho = 1, 0.6, 0.3$. +1. Explain the intuition. +``` + +```{solution-start} ta_ex4 +:class: dropdown +``` + +```{code-cell} ipython3 +def expected_utility(b, v_own, n, ρ): + "Expected utility of bidding b when rivals bid (n-1)v/(n-1+ρ)." + win_prob = np.minimum(1, b*(n - 1 + ρ)/(n - 1))**(n - 1) + return win_prob * np.maximum(v_own - b, 0)**ρ + +n_ex, v_own = 5, 0.8 +grid = np.linspace(0.001, v_own, 2001) + +print(f"{'ρ':>6}{'theory b*':>12}{'numerical':>12}") +for ρ in (1.0, 0.5, 0.2): + theory = (n_ex - 1)*v_own/(n_ex - 1 + ρ) + numerical = grid[int(np.argmax(expected_utility(grid, v_own, n_ex, ρ)))] + print(f"{ρ:>6}{theory:>12.4f}{numerical:>12.4f}") +``` + +```{code-cell} ipython3 +rng_ra = np.random.default_rng(42) +v_ra = np.sort(rng_ra.uniform(0, 1, (n_ex, 200_000)), axis=0) + +print(f"{'ρ':>6}{'FPSB revenue':>15}{'SPSB revenue':>15}") +for ρ in (1.0, 0.6, 0.3): + fpsb = ((n_ex - 1)/(n_ex - 1 + ρ)) * v_ra[-1, :] + print(f"{ρ:>6}{fpsb.mean():>15.4f}{v_ra[-2, :].mean():>15.4f}") +``` + +With risk neutrality ($\rho = 1$) the two auctions raise the same revenue, as the theorem says. + +With risk aversion ($\rho < 1$) the FPSB auction raises **more**. + +The intuition is that in a FPSB auction, shading one's bid is a gamble: it raises the surplus conditional on winning but lowers the probability of winning. + +A risk-averse bidder dislikes that gamble and so shades less, which transfers revenue to the seller. + +In a SPSB auction the winner's payment does not depend on her own bid, so risk aversion changes nothing: bidding one's valuation remains weakly dominant, and the seller's revenue is unaffected. + +```{solution-end} +``` + +## Further reading + +The second-price sealed-bid auction was proposed by {cite}`Vickrey_61`. + +For textbook treatments of the material in this lecture, see {cite}`Krishna2009` and {cite}`Milgrom2004`. + +{cite}`Klemperer1999` surveys the literature. + +The revenue equivalence theorem in the general form sketched above is due to {cite}`Myerson1981` and {cite}`RileySamuelson1981`. -1. Wikipedia for FPSB: https://en.wikipedia.org/wiki/First-price_sealed-bid_auction -2. Wikipedia for SPSB: https://en.wikipedia.org/wiki/Vickrey_auction -3. Chandra Chekuri's lecture note for algorithmic game theory: https://chekuri.cs.illinois.edu/teaching/spring2008/Lectures/scribed/Notes20.pdf -4. Tim Salmon. ECO 4400 Supplemental Handout: All About Auctions: https://s2.smu.edu/tsalmon/auctions.pdf -5. Auction Theory- Revenue Equivalence Theorem: https://michaellevet.wordpress.com/2015/07/06/auction-theory-revenue-equivalence-theorem/ -6. Order Statistics: https://online.stat.psu.edu/stat415/book/export/html/834 +Both auctions studied here are naturally described in terms of order statistics of bidders' valuations, a subject treated at length by {cite}`DavidNagaraja2003`.