Lucas Asset Pricing Model

A notebook by Christopher D. Carroll and Mateo Velásquez-Giraldo¶

Inspired by its Quantecon counterpart¶

This notebook presents simple computational tools to solve an instance of Lucas’s asset-pricing model for which there is no analytical solution: The case when the logarithm of the asset’s dividend follows an autoregressive process of order 1,

\ln d_{t+1} = \gamma + \alpha \ln d_t + \varepsilon_{t+1}, \qquad \varepsilon \sim \mathcal{N}(-\frac{\sigma^2}{2}, \sigma).

A presentation of this model can be found in Christopher D. Carroll’s lecture notes.

Those notes use the Bellman equation to derive a relationship between the price of the asset in the current period $t$ and the next period $t+1$ :

P_{t} = \overbrace{\left(\frac{1}{1+\vartheta}\right)} ^{\beta}\mathbb{E}_{t}\left[ \frac{u^{\prime}(d_{t+1})}{u^{\prime}(d_t)} (P_{t+1} + d_{t+1}) \right]

The equilibrium pricing equation is a relationship between the price and the dividend (a “pricing kernel”) $P^{*}(d)$ such that, if everyone believes that to be the pricing kernel, everyone’s Euler equation will be satisfied:

P^*(d_t) = \left(\frac{1}{1+\vartheta}\right)\mathbb{E}_{t}\left[ \frac{u^{\prime}(d_{t+1})}{u^{\prime}(d_t)} (P^*(d_{t+1}) + d_{t+1}) \right]

As noted in the handout, there are some special circumstances in which it is possible to solve for $P^{*}$ analytically:

Shock Process	Mean Restrictions	CRRA	Solution for Pricing Kernel $P^*(d)$
bounded, IID, $\mathbb{E}[d]=\bar{d}$	$0 < d < \infty$	1 (log)	$\vartheta^{-1}d$
lognormal, mean 1	$\mu=-\sigma^{2}/2$	$\rho$	$\vartheta^{-1}{d_t^\rho}~e^{\rho(\rho-1)\sigma^2/2}$
lognormal mean $e^{\gamma}=\mathbb{E}[d_{t+1}/d_{t}]$	${\mu~=-\sigma^{2}/2+\gamma}$	$\rho$	$\vartheta^{-1}d_t^{\rho}e^{\rho\gamma+\rho(\rho-1)\sigma^2/2}$

However, under most circumstances, the only way to obtain the pricing function $P^{*}$ is by solving for it numerically, as outlined below.

Finding the equilibrium pricing function.¶

We know that the equilibrium pricing function must satisfy the equation above. Let’s define an operator that allows us to evaluate whether any candidate pricing function satisfies this requirement.

Let $T$ be an operator which takes as argument a function and returns another function (these are usually called functionals or higher-order functions). For some function $f$ , denote with $T[f]$ the function that results from applying $T$ to $f$ . Then, for any real number $x$ , $T[f](x)$ will be the real number that one obtains when the function $T[f]$ is given $x$ as an input.

We define our particular operator as follows. For any function $g:\mathbb{R}\rightarrow\mathbb{R}$ , $T[g]$ is obtained as

\forall~d_t \in \mathbb{R},\,\,\,\, T[g](d_t) := \beta~\mathbb{E}_{t}\left[ \frac{u^{\prime}(d_{t+1})}{u^{\prime}(d_t)} (f(d_{t+1}) + d_{t+1}) \right].

We can use $T$ to re-express our pricing equation. If $P^*(\bullet)$ is our equilibrium pricing funtion, it must satisfy

\forall~d_t,\,\,\,\,P^*(d_t) = \beta\mathbb{E}_{t}\left[ \frac{u^{\prime}(d_{t+1})}{u^{\prime}(d_t)} (P^*(d_{t+1}) + d_{t+1}) \right] = T[P^*](d_t).

or, expressed differently,

P^* = T[P^*].

Our equilibrium pricing function is therefore a fixed point of the operator $T$ .

It turns out that $T$ is a contraction mapping. This is useful because it implies, through Banach’s fixed-point theorem, that:

$T$ has exactly one fixed point.
Starting from an arbitrary function $f$ , the sequence $\{T^n[f]\}_{n=1}^{\infty}$ converges to that fixed point.

For our purposes, this translates to:

Our equilibrium pricing function not only exists, but is unique.
We can get arbitrarily close to the equilibrium pricing function by making some initial guess $f$ and applying the operator $T$ to it repeatedly.

The code below creates a representation of our model and implements a solution routine to find $P^*$ . The main components of this routine are:

priceOnePeriod: this is operator $T$ from above. It takes a function $f$ , computes $\beta~\mathbb{E}_{t}\left[ \frac{u^{\prime}(d_{t+1})}{u^{\prime}(d_t)} (f(d_{t+1}) + d_{t+1}) \right]$ for a grid of $d_t$ values, and uses the result to construct a piecewise linear interpolator that approximates $T[f]$ .
solve: this is our iterative solution procedure. It generates an initial guess $f$ and applies priceOnePeriod to it iteratively. At each application, it constructs a measure of how much the candidate pricing function changed. Once changes between successive iterations are small enough, it declares that the solution has converged.

A computational representation of the problem and its solution.¶

Uninteresting setup:

# Setup
import numpy as np
import matplotlib.pyplot as plt
from copy import copy

from HARK.rewards import CRRAutilityP
from HARK.distributions import Normal, expected
from HARK.interpolation import LinearInterp, ConstantFunction

# A python class representing log-AR1 dividend processes.


class DivProcess:
    def __init__(self, α, σ, γ=0.0, nApprox=7):
        self.α = α
        self.σ = σ
        self.γ = γ
        self.nApprox = nApprox

        # Create a discrete approximation to the random shock
        self.ShkAppDstn = Normal(mu=-(σ**2) / 2, sigma=σ).discretize(N=nApprox)

    def getLogdGrid(self, n=100):
        """
        A method for creating a reasonable grid for log-dividends.
        """
        μ = self.γ - (self.σ**2) / 2
        uncond_sd = self.σ / np.sqrt(1 - self.α**2)
        uncond_mean = μ / (1 - self.α)
        logDGrid = np.linspace(-5 * uncond_sd, 5 * uncond_sd, n) + uncond_mean
        return logDGrid


# A class representing economies with Lucas trees.
class LucasEconomy:
    """
    A representation of an economy in which there are Lucas trees
    whose dividends' logarithm follows an AR1 process.
    """

    def __init__(self, CRRA, DiscFac, DivProcess):
        self.CRRA = CRRA
        self.DiscFac = DiscFac
        self.DivProcess = DivProcess
        self.uP = lambda c: CRRAutilityP(c, self.CRRA)

    def priceOnePeriod(self, Pfunc_next, logDGrid):
        # Create a function that, given current dividends
        # and the value of next period's shock, returns
        # the discounted value derived from the asset next period.
        def discounted_value(shock, log_d_now):
            # Find dividends
            d_now = np.exp(log_d_now)
            log_d_next = self.DivProcess.γ + self.DivProcess.α * log_d_now + shock
            d_next = np.exp(log_d_next)

            # Payoff and sdf
            payoff_next = Pfunc_next(log_d_next) + d_next
            SDF = self.DiscFac * self.uP(d_next / d_now)

            return SDF * payoff_next

        # The price at a given d_t is the expectation of the discounted value.
        # Compute it at every d in our grid. The expectation is taken over next
        # period's shocks
        prices_now = expected(
            dstn=self.DivProcess.ShkAppDstn, func=discounted_value, args=(logDGrid,)
        )

        # Create new interpolating price function
        Pfunc_now = LinearInterp(logDGrid, prices_now, lower_extrap=True)

        return Pfunc_now

    def solve(self, Pfunc_0=None, logDGrid=None, tol=1e-5, maxIter=500, disp=False):
        # Initialize the norm
        norm = tol + 1

        # Initialize Pfunc if initial guess is not provided
        if Pfunc_0 is None:
            Pfunc_0 = ConstantFunction(0.0)

        # Create a grid for log-dividends if one is not provided
        if logDGrid is None:
            logDGrid = self.DivProcess.getLogdGrid()

        # Initialize function and compute prices on the grid
        Pf_0 = copy(Pfunc_0)
        P_0 = Pf_0(logDGrid)

        it = 0
        while norm > tol and it < maxIter:
            # Apply the pricing equation
            Pf_next = self.priceOnePeriod(Pf_0, logDGrid)
            # Find new prices on the grid
            P_next = Pf_next(logDGrid)
            # Measure the change between price vectors
            norm = np.linalg.norm(P_0 - P_next)
            # Update price function and vector
            Pf_0 = Pf_next
            P_0 = P_next
            it = it + 1
            # Print iteration information
            if disp:
                print("Iter:" + str(it) + "   Norm = " + str(norm))

        if disp:
            if norm <= tol:
                print("Price function converged!")
            else:
                print("Maximum iterations exceeded!")

        self.EqlogPfun = Pf_0
        self.EqPfun = lambda d: self.EqlogPfun(np.log(d))

Creating and solving an example economy with AR1 dividends¶

An economy is fully specified by:

The dividend process for the assets (trees): we assume that $\ln d_{t+1} = \alpha \ln d_t + \varepsilon_{t+1}$ , $\varepsilon_{t+1}\sim\mathcal{N}(-\sigma^2/2,\sigma)$ . We must create a dividend process specifying $\alpha$ and $\sigma_{\varepsilon}$ .
The coefficient of relative risk aversion (CRRA).
The time-discount factor ( $\beta$ ).

# Create a log-AR1 process for dividends
DivProc = DivProcess(α=0.90, σ=0.1)

# Create an example economy
economy = LucasEconomy(CRRA=2, DiscFac=0.95, DivProcess=DivProc)

Once created, the economy can be ‘solved’, which means finding the equilibrium price kernel. The distribution of dividends at period $t+1$ depends on the value of dividends at $t$ , which also determines the resources agents have available to buy trees. Thus, $d_t$ is a state variable for the economy. The pricing function gives the price of trees that equates their demand and supply at every level of current dividends $d_t$ .

# Solve the economy, displaying the error term for each iteration
economy.solve(disp=True)

# After the economy is solved, we can use its Equilibrium price function
# to tell us the price if the dividend is 1
dvdnd = 1
print("P({}) = {:.6}".format(dvdnd, economy.EqPfun(dvdnd)))

Iter:1   Norm = 14.34310594086344
Iter:2   Norm = 14.614867281705973
Iter:3   Norm = 14.81462233119635
Iter:4   Norm = 14.939364249843898
Iter:5   Norm = 14.989619417063073
Iter:6   Norm = 14.968432855040843
Iter:7   Norm = 14.880572667620452
Iter:8   Norm = 14.731927474014178
Iter:9   Norm = 14.52902047558161
Iter:10   Norm = 14.278638854628227
Iter:11   Norm = 13.987555448796973
Iter:12   Norm = 13.662325970904508
Iter:13   Norm = 13.309147997993502
Iter:14   Norm = 12.933769279362002
Iter:15   Norm = 12.54143477811982
Iter:16   Norm = 12.136863389363548
Iter:17   Norm = 11.724246791512238
Iter:18   Norm = 11.307264248656217
Iter:19   Norm = 10.88910839449017
Iter:20   Norm = 10.472518082357068
Iter:21   Norm = 10.059815280329328
Iter:22   Norm = 9.652943734511805
Iter:23   Norm = 9.253507730992833
Iter:24   Norm = 8.862809773368175
Iter:25   Norm = 8.481886375493007
Iter:26   Norm = 8.111541464653394
Iter:27   Norm = 7.752377114027517
Iter:28   Norm = 7.40482148881162
Iter:29   Norm = 7.069154009563337
Iter:30   Norm = 6.745527819188153
Iter:31   Norm = 6.43398969488094
Iter:32   Norm = 6.134497580005358
Iter:33   Norm = 5.846935928796937
Iter:34   Norm = 5.571129063207857
Iter:35   Norm = 5.306852739535748
Iter:36   Norm = 5.053844115275684
Iter:37   Norm = 4.8118102958592255
Iter:38   Norm = 4.5804356280614105
Iter:39   Norm = 4.359387892926313
Iter:40   Norm = 4.148323536855
Iter:41   Norm = 3.946892065541599
Iter:42   Norm = 3.7547397120884534
Iter:43   Norm = 3.5715124781038123
Iter:44   Norm = 3.3968586350059407
Iter:45   Norm = 3.23043076218392
Iter:46   Norm = 3.071887389099951
Iter:47   Norm = 2.92089429983745
Iter:48   Norm = 2.777125550947055
Iter:49   Norm = 2.6402642466621504
Iter:50   Norm = 2.5100031095739235
Iter:51   Norm = 2.386044879600522
Iter:52   Norm = 2.2681025694850425
Iter:53   Norm = 2.155899601045683
Iter:54   Norm = 2.049169842908783
Iter:55   Norm = 1.9476575674272367
Iter:56   Norm = 1.8511173418645133
Iter:57   Norm = 1.7593138666580201
Iter:58   Norm = 1.672021771624408
Iter:59   Norm = 1.5890253792862712
Iter:60   Norm = 1.510118443058598
Iter:61   Norm = 1.4351038667929092
Iter:62   Norm = 1.3637934111186265
Iter:63   Norm = 1.2960073911146144
Iter:64   Norm = 1.2315743690705385
Iter:65   Norm = 1.1703308454400114
Iter:66   Norm = 1.1121209505264746
Iter:67   Norm = 1.0567961389685152
Iter:68   Norm = 1.004214888688016
Iter:69   Norm = 0.9542424056244561
Iter:70   Norm = 0.9067503352925078
Iter:71   Norm = 0.8616164819572083
Iter:72   Norm = 0.818724536020655
Iter:73   Norm = 0.777963810043275
Iter:74   Norm = 0.7392289836834438
Iter:75   Norm = 0.7024198577214298
Iter:76   Norm = 0.6674411172382658
Iter:77   Norm = 0.634202103941166
Iter:78   Norm = 0.6026165975626904
Iter:79   Norm = 0.5726026062095858
Iter:80   Norm = 0.5440821654963686
Iter:81   Norm = 0.5169811462657242
Iter:82   Norm = 0.49122907067358845
Iter:83   Norm = 0.4667589363979194
Iter:84   Norm = 0.44350704871700236
Iter:85   Norm = 0.421412860193699
Iter:86   Norm = 0.4004188176968008
Iter:87   Norm = 0.3804702164882355
Iter:88   Norm = 0.3615150611045491
Iter:89   Norm = 0.3435039327628806
Iter:90   Norm = 0.32638986302583317
Iter:91   Norm = 0.31012821346309866
Iter:92   Norm = 0.2946765610546749
Iter:93   Norm = 0.27999458908628433
Iter:94   Norm = 0.2660439832949772
Iter:95   Norm = 0.2527883330301617
Iter:96   Norm = 0.24019303720337631
Iter:97   Norm = 0.2282252148079118
Iter:98   Norm = 0.21685361979727616
Iter:99   Norm = 0.20604856012013625
Iter:100   Norm = 0.19578182071641212
Iter:101   Norm = 0.18602659028839066
Iter:102   Norm = 0.1767573916672117
Iter:103   Norm = 0.16795001560349337
Iter:104   Norm = 0.15958145781814148
Iter:105   Norm = 0.1516298591567478
Iter:106   Norm = 0.14407444869702674
Iter:107   Norm = 0.13689548966714943
Iter:108   Norm = 0.1300742280379809
Iter:109   Norm = 0.1235928436590115
Iter:110   Norm = 0.1174344038138782
Iter:111   Norm = 0.11158281907722094
Iter:112   Norm = 0.10602280136006649
Iter:113   Norm = 0.10073982403589372
Iter:114   Norm = 0.09572008404579596
Iter:115   Norm = 0.09095046588484632
Iter:116   Norm = 0.08641850737699282
Iter:117   Norm = 0.08211236714998044
Iter:118   Norm = 0.07802079372731051
Iter:119   Norm = 0.07413309615584784
Iter:120   Norm = 0.07043911609430371
Iter:121   Norm = 0.06692920128960075
Iter:122   Norm = 0.06359418037268635
Iter:123   Norm = 0.06042533890809264
Iter:124   Norm = 0.057414396635407663
Iter:125   Norm = 0.054553485843217786
Iter:126   Norm = 0.05183513081936258
Iter:127   Norm = 0.049252228324352595
Iter:128   Norm = 0.04679802903643489
Iter:129   Norm = 0.04446611992096433
Iter:130   Norm = 0.04225040747722925
Iter:131   Norm = 0.040145101819746856
Iter:132   Norm = 0.03814470155217281
Iter:133   Norm = 0.03624397939444744
Iter:134   Norm = 0.03443796852609161
Iter:135   Norm = 0.03272194960916258
Iter:136   Norm = 0.031091438458342148
Iter:137   Norm = 0.029542174324315832
Iter:138   Norm = 0.028070108761445414
Iter:139   Norm = 0.02667139504984569
Iter:140   Norm = 0.025342378144178373
Iter:141   Norm = 0.024079585123435457
Iter:142   Norm = 0.022879716116425367
Iter:143   Norm = 0.021739635679333864
Iter:144   Norm = 0.020656364602714354
Iter:145   Norm = 0.019627072127070527
Iter:146   Norm = 0.018649068545766374
Iter:147   Norm = 0.01771979817680735
Iter:148   Norm = 0.016836832684830137
Iter:149   Norm = 0.015997864735649348
Iter:150   Norm = 0.015200701966844189
Iter:151   Norm = 0.014443261259211396
Iter:152   Norm = 0.013723563293191268
Iter:153   Norm = 0.013039727376824886
Iter:154   Norm = 0.012389966531120225
Iter:155   Norm = 0.011772582820757146
Iter:156   Norm = 0.011185962916935679
Iter:157   Norm = 0.01062857388152953
Iter:158   Norm = 0.010098959161571077
Iter:159   Norm = 0.009595734782801998
Iter:160   Norm = 0.00911758573331583
Iter:161   Norm = 0.008663262527145959
Iter:162   Norm = 0.008231577939258631
Iter:163   Norm = 0.007821403903143986
Iter:164   Norm = 0.007431668562953552
Iter:165   Norm = 0.007061353472598732
Iter:166   Norm = 0.006709490934464021
Iter:167   Norm = 0.0063751614704972495
Iter:168   Norm = 0.006057491419614152
Iter:169   Norm = 0.005755650654652022
Iter:170   Norm = 0.005468850413064407
Iter:171   Norm = 0.005196341235730253
Iter:172   Norm = 0.004937411008699256
Iter:173   Norm = 0.004691383101932497
Iter:174   Norm = 0.004457614601632301
Iter:175   Norm = 0.004235494629770892
Iter:176   Norm = 0.004024442748025629
Iter:177   Norm = 0.003823907440919708
Iter:178   Norm = 0.0036333646745771654
Iter:179   Norm = 0.003452316527362579
Iter:180   Norm = 0.0032802898888329274
Iter:181   Norm = 0.0031168352231918865
Iter:182   Norm = 0.002961525394785437
Iter:183   Norm = 0.0028139545518366246
Iter:184   Norm = 0.00267373706579672
Iter:185   Norm = 0.0025405065238508634
Iter:186   Norm = 0.0024139147712710584
Iter:187   Norm = 0.0022936310014923895
Iter:188   Norm = 0.0021793408920568543
Iter:189   Norm = 0.0020707457828244703
Iter:190   Norm = 0.0019675618956923836
Iter:191   Norm = 0.0018695195930966792
Iter:192   Norm = 0.0017763626732945814
Iter:193   Norm = 0.0016878477009181062
Iter:194   Norm = 0.001603743370814542
Iter:195   Norm = 0.001523829903624448
Iter:196   Norm = 0.0014478984714075974
Iter:197   Norm = 0.001375750651995658
Iter:198   Norm = 0.0013071979105175545
Iter:199   Norm = 0.001242061106635624
Iter:200   Norm = 0.0011801700264542725
Iter:201   Norm = 0.0011213629376277934
Iter:202   Norm = 0.0010654861670085145
Iter:203   Norm = 0.0010123936987695278
Iter:204   Norm = 0.0009619467929445259
Iter:205   Norm = 0.0009140136229716336
Iter:206   Norm = 0.0008684689310220655
Iter:207   Norm = 0.0008251937008250442
Iter:208   Norm = 0.0007840748465884144
Iter:209   Norm = 0.0007450049175648438
Iter:210   Norm = 0.0007078818171773037
Iter:211   Norm = 0.0006726085362432564
Iter:212   Norm = 0.0006390928994664652
Iter:213   Norm = 0.0006072473245952711
Iter:214   Norm = 0.0005769885936040177
Iter:215   Norm = 0.0005482376350572404
Iter:216   Norm = 0.0005209193177180375
Iter:217   Norm = 0.0004949622539685552
Iter:218   Norm = 0.00047029861347630916
Iter:219   Norm = 0.0004468639458267492
Iter:220   Norm = 0.0004245970121032102
Iter:221   Norm = 0.00040343962488112255
Iter:222   Norm = 0.00038333649621035825
Iter:223   Norm = 0.00036423509313689677
Iter:224   Norm = 0.00034608550026885444
Iter:225   Norm = 0.00032884028958809217
Iter:226   Norm = 0.0003124543962759582
Iter:227   Norm = 0.0002968850011565627
Iter:228   Norm = 0.00028209141866892884
Iter:229   Norm = 0.0002680349905576475
Iter:230   Norm = 0.0002546789849357291
Iter:231   Norm = 0.00024198850023253644
Iter:232   Norm = 0.0002299303738882025
Iter:233   Norm = 0.00021847309604241106
Iter:234   Norm = 0.00020758672670371484
Iter:235   Norm = 0.000197242817963524
Iter:236   Norm = 0.00018741433928826714
Iter:237   Norm = 0.00017807560717920262
Iter:238   Norm = 0.00016920221791301237
Iter:239   Norm = 0.00016077098373458665
Iter:240   Norm = 0.00015275987235465702
Iter:241   Norm = 0.00014514794933610434
Iter:242   Norm = 0.000137915323380576
Iter:243   Norm = 0.00013104309440237102
Iter:244   Norm = 0.00012451330400711154
Iter:245   Norm = 0.00011830888873856297
Iter:246   Norm = 0.00011241363536902475
Iter:247   Norm = 0.0001068121385326778
Iter:248   Norm = 0.00010148976060521523
Iter:249   Norm = 9.643259320029383e-05
Iter:250   Norm = 9.162742108566186e-05
Iter:251   Norm = 8.706168754619847e-05
Iter:252   Norm = 8.272346145365097e-05
Iter:253   Norm = 7.860140630130716e-05
Iter:254   Norm = 7.468475039819174e-05
Iter:255   Norm = 7.096325889746349e-05
Iter:256   Norm = 6.742720682601107e-05
Iter:257   Norm = 6.406735387874806e-05
Iter:258   Norm = 6.087492015769238e-05
Iter:259   Norm = 5.784156329083827e-05
Iter:260   Norm = 5.495935660985967e-05
Iter:261   Norm = 5.222076833999482e-05
Iter:262   Norm = 4.9618642107046365e-05
Iter:263   Norm = 4.7146178127042934e-05
Iter:264   Norm = 4.4796915380576354e-05
Iter:265   Norm = 4.256471482312034e-05
Iter:266   Norm = 4.044374335343422e-05
Iter:267   Norm = 3.8428458452395206e-05
Iter:268   Norm = 3.65135939085278e-05
Iter:269   Norm = 3.469414577592586e-05
Iter:270   Norm = 3.296535953101107e-05
Iter:271   Norm = 3.132271751942793e-05
Iter:272   Norm = 2.9761927333477716e-05
Iter:273   Norm = 2.827891023551752e-05
Iter:274   Norm = 2.686979091958857e-05
Iter:275   Norm = 2.5530887075460387e-05
Iter:276   Norm = 2.4258699903291916e-05
Iter:277   Norm = 2.304990500900857e-05
Iter:278   Norm = 2.1901343559741297e-05
Iter:279   Norm = 2.0810014136204545e-05
Iter:280   Norm = 1.9773064945727455e-05
Iter:281   Norm = 1.8787786251643173e-05
Iter:282   Norm = 1.7851603329965908e-05
Iter:283   Norm = 1.696206980654742e-05
Iter:284   Norm = 1.611686111768278e-05
Iter:285   Norm = 1.5313768645419402e-05
Iter:286   Norm = 1.4550693695290708e-05
Iter:287   Norm = 1.3825642320721789e-05
Iter:288   Norm = 1.313671972337381e-05
Iter:289   Norm = 1.2482125690853726e-05
Iter:290   Norm = 1.1860149638043228e-05
Iter:291   Norm = 1.126916627003413e-05
Iter:292   Norm = 1.0707631213278965e-05
Iter:293   Norm = 1.0174077062486109e-05
Iter:294   Norm = 9.667109531709195e-06
Price function converged!
P(1) = 20.1571

The effect of risk aversion.¶

The notes discuss the surprising implication that an increase in the coefficient of relative risk aversion $\rho$ leads to higher prices for the risky trees! This is demonstrated below.

# Create two economies with different risk aversion
Disc = 0.95
LowCRRAEcon = LucasEconomy(CRRA=2, DiscFac=Disc, DivProcess=DivProc)
HighCRRAEcon = LucasEconomy(CRRA=4, DiscFac=Disc, DivProcess=DivProc)

# Solve both
LowCRRAEcon.solve()
HighCRRAEcon.solve()

# Plot the pricing functions for both
dGrid = np.linspace(0.5, 2.5, 30)
plt.figure()
plt.plot(dGrid, LowCRRAEcon.EqPfun(dGrid), label="Low CRRA")
plt.plot(dGrid, HighCRRAEcon.EqPfun(dGrid), label="High CRRA")
plt.legend()
plt.xlabel("$d_t$")
plt.ylabel("$P_t$")

Testing our analytical solutions¶

Case 1: Log Utility¶

The lecture notes show that with logarithmic utility (a CRRA of 1), the pricing kernel has a closed form expression:

P^*(d_t) = \frac{d_t}{\vartheta}

(1)

We now compare our numerical solution with this analytical expression.

# Create an economy with log utility and the same dividend process from before
logUtilEcon = LucasEconomy(CRRA=1, DiscFac=Disc, DivProcess=DivProc)
# Solve it
logUtilEcon.solve()

# Generate a function with our analytical solution
theta = 1 / Disc - 1


def aSol(d):
    return d / theta


# Get a grid for d over which to compare them
dGrid = np.exp(DivProc.getLogdGrid())

# Plot both
plt.figure()
plt.plot(dGrid, aSol(dGrid), "*", label="Analytical solution")
plt.plot(dGrid, logUtilEcon.EqPfun(dGrid), label="Numerical solution")
plt.legend()
plt.xlabel("$d_t$")
plt.ylabel("$P^*(d_t)$")

Case 2: I.I.D dividends¶

The notes also show that, if $\ln d_{t+n}\sim \mathcal{N}(-\sigma^2/2, \sigma^2)$ for all $n$ , the pricing kernel is exactly

P^*(d_t) = d_t^\rho\times e^{\rho(\rho-1)\sigma^2/2}\overbrace{\frac{\beta}{1-\beta}}^{=\vartheta}

We now our numerical solution for this case.

# Create an i.i.d. dividend process
σ = 0.1
iidDivs = DivProcess(α=0.0, σ=σ)

# And an economy that embeds it
CRRA = 2
Disc = 0.9

iidEcon = LucasEconomy(CRRA=CRRA, DiscFac=Disc, DivProcess=iidDivs)
iidEcon.solve()

# Generate a function with our analytical solution
dTil = np.exp((σ**2) / 2 * CRRA * (CRRA - 1))


def aSolIID(d):
    return d**CRRA * dTil * Disc / (1 - Disc)


# Get a grid for d over which to compare them
dGrid = np.exp(iidDivs.getLogdGrid())

# Plot both
plt.figure()
plt.plot(dGrid, aSolIID(dGrid), "*", label="Analytical solution")
plt.plot(dGrid, iidEcon.EqPfun(dGrid), label="Numerical solution")
plt.legend()
plt.xlabel("$d_t$")
plt.ylabel("$P^*(d_t)$")
plt.show()

Case 3: Dividends that are a geometric random walk with drift¶

The notes also show that if the dividend process is

\ln d_{t+1} = \gamma + \ln d_t + \varepsilon_{t+1}, \qquad \varepsilon \sim \mathcal{N}(-\frac{\sigma^2}{2}, \sigma).

so that $E_t[d_{t+1}/d_t] = e^\gamma$ , then we have

P^*(d_t) = d_t^\rho\times e^{(\rho-1)\left(\rho\sigma^2/2 - \gamma\right)}\overbrace{\frac{\beta}{1-\beta}}^{\vartheta}

which, when $\rho=1$ , reduces (as it should) to

\frac{P^*(d_t)}{d_t} = \vartheta

CRRA = 2
Disc = 0.9
σ = 0.1
γ = 0.3

# Create a random walk dividend process
# (it turns out that the whole model can be normalized by d_t, and
# in normalized, terms the dividend proces becomes iid again)
rw_divs = DivProcess(γ=γ, α=0, σ=σ)

# And an economy that embeds it
rw_econ = LucasEconomy(CRRA=CRRA, DiscFac=Disc, DivProcess=rw_divs)
rw_econ.solve()

# Generate a function with our analytical solution
a_sol_factor = np.exp((CRRA - 1) * (CRRA * σ**2 / 2 - γ))


def a_sol_rw(d):
    return d**CRRA * a_sol_factor * Disc / (1 - Disc)


# Get a grid for d over which to compare them
dGrid = np.exp(rw_divs.getLogdGrid())

# Plot both
plt.figure()
plt.plot(dGrid, a_sol_rw(dGrid), "*", label="Analytical solution")
plt.plot(dGrid, rw_econ.EqPfun(dGrid), label="Numerical solution")
plt.legend()
plt.xlabel("$d_t$")
plt.ylabel("$P^*(d_t)$")
plt.show()

Testing our approximation of the dividend process¶

Hidden in the solution method implemented above is the fact that, in order to make expectations easy to compute, we discretize the random shock $\varepsilon_t$ , which is to say, we create a discrete variable $\tilde{\varepsilon}$ that approximates the behavior of $\varepsilon_t$ . This is done using a Gauss-Hermite quadrature.

A parameter for the numerical solution is the number of different values that we allow our discrete approximation $\tilde{\varepsilon}$ to take, $n^{\#}$ . We would expect a higher n^# to improve our solution, as the discrete approximation of $\varepsilon_t$ improves. We test this below.

# Increase CRRA to make the effect of uncertainty more evident.
CRRA = 10
Disc = 0.9
σ = 0.1
ns = [1, 2, 10]

dTil = np.exp((σ**2) / 2 * CRRA * (CRRA - 1.0))


def aSolIID(d):
    return d**CRRA * dTil * Disc / (1 - Disc)


dGrid = np.exp(iidDivs.getLogdGrid())

plt.figure()
for n in ns:
    iidDivs = DivProcess(α=0.0, σ=σ, nApprox=n)
    iidEcon = LucasEconomy(CRRA=CRRA, DiscFac=Disc, DivProcess=iidDivs)
    iidEcon.solve()
    plt.plot(dGrid, iidEcon.EqPfun(dGrid), label="Num.Sol. $n^\#$ = {}".format(n))

# Plot both
plt.plot(dGrid, aSolIID(dGrid), "*", label="Analytical solution")
plt.legend()
plt.xlabel("$d_t$")
plt.ylabel("$P^*(d_t)$")
plt.show()