<jupyter_start><jupyter_text>Part 3.1 Python Good Code, OOP I and OOP II  An example of Bad CodeLet’s have a look at some poorly written codeThe job of the code is to generate and plot time series of the simplified Solow model$k_{t+1} = s k_t^{\alpha} + (1 - \delta) k_t,\quad t = 0, 1, 2, \ldots$Here* $kt$ is capital at time $t$ and* $s,α,δ$ are parameters (savings, a productivity parameter and depreciation)For each parameterization, the code1. sets $k0$=12. iterates using (1) to produce a sequence $k_0,k_1,k_2…,k_T$3. plots the sequenceThe plots will be grouped into three subfiguresIn each subfigure, two parameters are held fixed while another varies<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

# Allocate memory for time series
k = np.empty(50)

fig, axes = plt.subplots(3, 1, figsize=(12, 15))

# Trajectories with different α
δ = 0.1
s = 0.4
α = (0.25, 0.33, 0.45)

for j in range(3):
    k[0] = 1
    for t in range(49):
        k[t+1] = s * k[t]**α[j] + (1 - δ) * k[t]
    axes[0].plot(k, 'o-', label=rf"$\alpha = {α[j]},\; s = {s},\; \delta={δ}$")

axes[0].grid(lw=0.2)
axes[0].set_ylim(0, 18)
axes[0].set_xlabel('time')
axes[0].set_ylabel('capital')
axes[0].legend(loc='upper left', frameon=True, fontsize=14)

# Trajectories with different s
δ = 0.1
α = 0.33
s = (0.3, 0.4, 0.5)

for j in range(3):
    k[0] = 1
    for t in range(49):
        k[t+1] = s[j] * k[t]**α + (1 - δ) * k[t]
    axes[1].plot(k, 'o-', label=rf"$\alpha = {α},\; s = {s},\; \delta={δ}$")

axes[1].grid(lw=0.2)
axes[1].set_xlabel('time')
axes[1].set_ylabel('capital')
axes[1].set_ylim(0, 18)
axes[1].legend(loc='upper left', frameon=True, fontsize=14)

# Trajectories with different δ
δ = (0.05, 0.1, 0.15)
α = 0.33
s = 0.4

for j in range(3):
    k[0] = 1
    for t in range(49):
        k[t+1] = s * k[t]**α + (1 - δ[j]) * k[t]
    axes[2].plot(k, 'o-', label=rf"$\alpha = {α},\; s = {s},\; \delta={δ[j]}$")

axes[2].set_ylim(0, 18)
axes[2].set_xlabel('time')
axes[2].set_ylabel('capital')
axes[2].grid(lw=0.2)
axes[2].legend(loc='upper left', frameon=True, fontsize=14)

plt.show()<jupyter_output><empty_output><jupyter_text>The code follows **PEP8**, but it's still bad code. Why?* It uses Magic Numbers, i.e. numeric literals in the code body instead of named constants. This hampers reusability, and causes readability and maintainability problems.* Blocks of logic are repeated, with minor changes, instead of taking advantage of the programming language, and automating it.* Makes liberal use of Global Variables, at the cost of code-determinism, and prevents the code from being efficiently/safely JIT-ted.<jupyter_code># A better example

from itertools import product

def plot_path(ax, αs, s_vals, δs, series_length=50):
    """
    Add a time series plot to the axes ax for all given parameters.
    """
    k = np.empty(series_length)

    for (α, s, δ) in product(αs, s_vals, δs):
        k[0] = 1
        for t in range(series_length-1):
            k[t+1] = s * k[t]**α + (1 - δ) * k[t]
        ax.plot(k, 'o-', label=rf"$\alpha = {α},\; s = {s},\; \delta = {δ}$")

    ax.grid(lw=0.2)
    ax.set_xlabel('time')
    ax.set_ylabel('capital')
    ax.set_ylim(0, 18)
    ax.legend(loc='upper left', frameon=True, fontsize=14)

fig, axes = plt.subplots(3, 1, figsize=(12, 15))

# Parameters (αs, s_vals, δs)
set_one = ([0.25, 0.33, 0.45], [0.4], [0.1])
set_two = ([0.33], [0.3, 0.4, 0.5], [0.1])
set_three = ([0.33], [0.4], [0.05, 0.1, 0.15])

for (ax, params) in zip(axes, (set_one, set_two, set_three)):
    αs, s_vals, δs = params
    plot_path(ax, αs, s_vals, δs)

plt.show()<jupyter_output><empty_output><jupyter_text>Using ClassesClasses are blueprints that help you build objects according to your own specifications.In the OOP paradigm, data and functions are bundled together into “objects”An example is a Python list, which not only stores data, but also knows how to sort itself, etc.```Pythonx = [1, 5, 4]x.sort()x``````None[1, 4, 5]````sort` is a function that is “part of” the list object — and hence called a *method*If we want to make our own types of objects we need to use class definitionsA class definition is a blueprint for a particular class of objects (e.g., lists, strings or complex numbers)It describes* What kind of data the class stores* What methods it has for acting on these dataAn object or instance is a realization of the class, created from the blueprint* Each instance has its own unique data* Methods set out in the class definition act on this (and other) dataIn Python, the data and methods of an object are collectively referred to as attributesAttributes are accessed via “dotted attribute notation”* `object_name.data`* `object_name.method_name()` Why OOP and how is it useful ?OOP is useful for the same reason that **abstraction** is useful: *for recognizing and exploiting common structures*For example,* a Markov chain consists of a set of states and a collection of transition probabilities for moving across states* a general equilibrium theory consists of a commodity space, preferences, technologies, and an equilibrium definition* a game consists of a list of players, lists of actions available to each player, player payoffs as functions of all players’ actions, and a timing protocolThese are all abstractions that collect together “objects” of the same “type”Recognizing common structure allows us to employ common toolsIn economic theory, this might be a proposition that applies to all games of a certain typeIn Python, this might be a method that’s useful for all Markov chains (e.g., `simulate`)When we use OOP, the `simulate` method is conveniently bundled together with the Markov chain object<jupyter_code># An example:

class Consumer:
    
    def __init(self:object, w:float)->None:
        self.wealth = w
        
    def earn(self, y:float)->None:
        self.wealth += y
        
    def spend(self, x:float)->None:
        
        if (self.wealth-x) < 0:
            print("Insufficient Funds")
        else:
            self.wealth -= x<jupyter_output><empty_output><jupyter_text>Example: The Solow Growth ModelFor our next example, let’s write a simple class to implement the Solow growth modelThe Solow growth model is a neoclassical growth model where the amount of capital stock per capita $kt$ evolves according to the rule$k_{t+1} = \frac{s z k_t^{\alpha} + (1 - \delta) k_t}{1 + n}$Here* $s$ is an exogenously given savings rate* $z$ is a productivity parameter* $α$ is capital’s share of income* $n$ is the population growth rate* $δ$ is the depreciation rateThe steady state of the model is the $k$ that solves when $kt+1=kt=k$Here’s a class that implements this modelSome points of interest in the code are* An instance maintains a record of its current capital stock in the variable self.k* The h method implements the right hand side of (1)* The update method uses h to update capital as per (1)    * Notice how inside update the reference to the local method `h` is `self.h`The methods `steady_state` and `generate_sequence` are fairly self explanatory<jupyter_code># Another Example

class Solow:
    
    def __init__(self,
                 n:float=0.05,  # population growth rate
                 s:float=0.25,  # savings rate
                 δ:float=0.1,  # depreciation rate
                 α:float=0.3,  # share of labor
                 z:float=2.0,  # productivity
                 k:float=1.0  # current capital stock
                )->None:
        
        self.nx = n
        self.sx = s
        self.dx = δ
        self.ax = α
        self.zx = z
        self.kx = k
        
    def h_func(self)->float:
        
        nx,sx,dx,ax,zx = self.nx,self.sx,self.dx,self.ax,self.zx
        
        return (sx * zx * (self.kx**ax) + (1 - dx) * self.kx) / (1 + nx)
    
    def update(self)->None:
        
        self.k = self.h_func()
    
    def steady_state(self)->float:
        
        nx,sx,dx,ax,zx = self.nx,self.sx,self.dx,self.ax,self.zx
        
        return ((sx * zx) / (nx + dx))**(1 / (1 - ax))
    
    def gen_seq(self, t:int)->list:
        
        path = []
        for i in range(t):
            path.append(self.kx)
            self.update()
        return path
import matplotlib.pyplot as plt

s1 = Solow()
s2 = Solow(k=8.0)

T = 60
fig, ax = plt.subplots(figsize=(9, 6))

# Plot the common steady state value of capital
ax.plot([s1.steady_state()]*T, 'k-', label='steady state')

# Plot time series for each economy
for s in s1, s2:
    lb = f'capital series from initial state {s.kx}'
    ax.plot(s.gen_seq(T), 'o-', lw=2, alpha=0.6, label=lb)

ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>EXample: A MarketNext let’s write a class for a simple one good market where agents are price takersThe market consists of the following objects:* A linear demand curve $Q= a_d − b_d p$* A linear supply curve $Q= a_z + b_z(p−t)$Here $p$ is price paid by the consumer, $Q$ is quantity, and $t$ is a per unit tax    * Other symbols are demand and supply parametersThe class provides methods to compute various values of interest, including competitive equlibrium price and quantity, tax revenue raised, consumer surplus and producer surplus<jupyter_code>from scipy.integrate import quad

class Market:

    def __init__(self, ad, bd, az, bz, tax):
        """
        Set up market parameters.  All parameters are scalars.  See
        https://lectures.quantecon.org/py/python_oop.html for interpretation.

        """
        self.ad, self.bd, self.az, self.bz, self.tax = ad, bd, az, bz, tax
        if ad < az:
            raise ValueError('Insufficient demand.')

    def price(self):
        "Return equilibrium price"
        return  (self.ad - self.az + self.bz * self.tax) / (self.bd + self.bz)

    def quantity(self):
        "Compute equilibrium quantity"
        return  self.ad - self.bd * self.price()

    def consumer_surp(self):
        "Compute consumer surplus"
        # == Compute area under inverse demand function == #
        integrand = lambda x: (self.ad / self.bd) - (1 / self.bd) * x
        area, error = quad(integrand, 0, self.quantity())
        return area - self.price() * self.quantity()

    def producer_surp(self):
        "Compute producer surplus"
        #  == Compute area above inverse supply curve, excluding tax == #
        integrand = lambda x: -(self.az / self.bz) + (1 / self.bz) * x
        area, error = quad(integrand, 0, self.quantity())
        return (self.price() - self.tax) * self.quantity() - area

    def taxrev(self):
        "Compute tax revenue"
        return self.tax * self.quantity()

    def inverse_demand(self, x):
        "Compute inverse demand"
        return self.ad / self.bd - (1 / self.bd)* x

    def inverse_supply(self, x):
        "Compute inverse supply curve"
        return -(self.az / self.bz) + (1 / self.bz) * x + self.tax

    def inverse_supply_no_tax(self, x):
        "Compute inverse supply curve without tax"
        return -(self.az / self.bz) + (1 / self.bz) * x

#sample of usage
baseline_params = 15, .5, -2, .5, 3
m = Market(*baseline_params)
print("equilibrium price = ", m.price())

print("consumer surplus = ", m.consumer_surp())<jupyter_output>equilibrium price =  18.5
consumer surplus =  33.0625<jupyter_text>Here’s a short program that uses this class to plot an inverse demand curve together with inverse supply curves with and without taxes<jupyter_code>import numpy as np

# Baseline ad, bd, az, bz, tax
baseline_params = 15, .5, -2, .5, 3
m = Market(*baseline_params)

q_max = m.quantity() * 2
q_grid = np.linspace(0.0, q_max, 100)
pd = m.inverse_demand(q_grid)
ps = m.inverse_supply(q_grid)
psno = m.inverse_supply_no_tax(q_grid)

fig, ax = plt.subplots()
ax.plot(q_grid, pd, lw=2, alpha=0.6, label='demand')
ax.plot(q_grid, ps, lw=2, alpha=0.6, label='supply')
ax.plot(q_grid, psno, '--k', lw=2, alpha=0.6, label='supply without tax')
ax.set_xlabel('quantity', fontsize=14)
ax.set_xlim(0, q_max)
ax.set_ylabel('price', fontsize=14)
ax.legend(loc='lower right', frameon=False, fontsize=14)
plt.show()<jupyter_output><empty_output><jupyter_text>The next program provides a function that* takes an instance of `Market` as a parameter* computes dead weight loss from the imposition of the tax<jupyter_code>def deadw(m):
    "Computes deadweight loss for market m."
    # == Create analogous market with no tax == #
    m_no_tax = Market(m.ad, m.bd, m.az, m.bz, 0)
    # == Compare surplus, return difference == #
    surp1 = m_no_tax.consumer_surp() + m_no_tax.producer_surp()
    surp2 = m.consumer_surp() + m.producer_surp() + m.taxrev()
    return surp1 - surp2
# Here’s an example of usage

baseline_params = 15, .5, -2, .5, 3
m = Market(*baseline_params)
deadw(m)  # Show deadweight loss<jupyter_output><empty_output><jupyter_text>Example: ChaosLet’s look at one more example, related to chaotic dynamics in nonlinear systemsOne simple transition rule that can generate complex dynamics is the logistic map$ x_{t+1}=rx_t(1−x_t),x_0∈[0,1],r∈[0,4] $Let’s write a class for generating time series from this modelHere’s one implementation<jupyter_code>class Chaos:
  """
  Models the dynamical system with :math:`x_{t+1} = r x_t (1 - x_t)`
  """
  def __init__(self, x0, r):
      """
      Initialize with state x0 and parameter r
      """
      self.x, self.r = x0, r

  def update(self):
      "Apply the map to update state."
      self.x =  self.r * self.x *(1 - self.x)

  def generate_sequence(self, n):
      "Generate and return a sequence of length n."
      path = []
      for i in range(n):
          path.append(self.x)
          self.update()
      return path
# Usage
ch = Chaos(0.1, 4.0)     # x0 = 0.1 and r = 0.4
ch.generate_sequence(5)  # First 5 iterates



ch = Chaos(0.1, 4.0)
ts_length = 250

fig, ax = plt.subplots()
ax.set_xlabel('$t$', fontsize=14)
ax.set_ylabel('$x_t$', fontsize=14)
x = ch.generate_sequence(ts_length)
ax.plot(range(ts_length), x, 'bo-', alpha=0.5, lw=2, label='$x_t$')
plt.show()<jupyter_output><empty_output><jupyter_text>The next piece of code provides a bifurcation diagram<jupyter_code>fig, ax = plt.subplots()
ch = Chaos(0.1, 4)
r = 2.5
while r < 4:
    ch.r = r
    t = ch.generate_sequence(1000)[950:]
    ax.plot([r] * len(t), t, 'b.', ms=0.6)
    r = r + 0.005

ax.set_xlabel('$r$', fontsize=16)
plt.show()<jupyter_output><empty_output><jupyter_text>On the horizontal axis is the parameter $r$The vertical axis is the state space $[0,1]$For each $r$ we compute a long time series and then plot the tail (the last 50 points)The tail of the sequence shows us where the trajectory concentrates after settling down to some kind of steady state, if a steady state exists.Whether it settles down, and the character of the steady state to which it does settle down, depend on the value of $r$For $r$ between about 2.5 and 3, the time series settles into a single fixed point plotted on the vertical axisFor $r$ between about 3 and 3.45, the time series settles down to oscillating between the two values plotted on the vertical axisFor $r$ a little bit higher than 3.45, the time series settles down to oscillating among the four values plotted on the vertical axisNotice that there is no value of $r$ that leads to a steady state oscillating among three values  Exercises Exercise 1The empirical cumulative distribution function (ecdf) corresponding to a sample $\{X_i\}_{i=1}^n$ is defined as * $F_n(x) := \frac{1}{n}  \sum_{i=1}^n \mathbf{1}\{X_i \leq x\}  \qquad (x \in \mathbb{R})$Here $1{Xi≤x}$ is an indicator function (one if $Xi≤x$ and zero otherwise) and hence $Fn(x)$ is the fraction of the sample that falls below $x$The Glivenko–Cantelli Theorem states that, provided that the sample is iid, the ecdf $F_n$ converges to the true distribution function $F$**Implement $Fn$** as a class called ECDF, where* A given sample $\{X_i\}_{i=1}^n$ are the instance data, stored as self.observations* The class implements a `__call__` method that returns $Fn(x)$for any $x$Your code should work as follows (modulo randomness)```Pythonfrom random import uniformsamples = [uniform(0, 1) for i in range(10)]F = ECDF(samples)F(0.5)   Evaluate ecdf at x = 0.5``````None0.29``````PythonF.observations = [uniform(0, 1) for i in range(1000)]F(0.5)``````None0.479```Aim for clarity, not efficiency<jupyter_code>class ECDF:
    
    #use no imports
    
    def __init__(self, observations:list)->None:
        
        self.obs = observations

    def __call__(self, x:float)->float:
        
        counter = len([1 for i in self.obs if i < x])
        return counter / len(self.obs)
        
from random import uniform

samples = [uniform(0, 1) for i in range(10)]
F = ECDF(samples)

print(F(0.5))  # Evaluate ecdf at x = 0.5

F.obs = [uniform(0, 1) for i in range(1000)]

print(F(0.5))<jupyter_output>0.4
0.473<jupyter_text>Exercise 2In an earlier exercise, you wrote a function for evaluating polynomialsThis exercise is an extension, where the task is to build a simple class called `Polynomial` for representing and manipulating polynomial functions such as$p(x) = a_0 + a_1 x + a_2 x^2 + \cdots a_N x^N = \sum_{n=0}^N a_n x^n    \qquad (x \in \mathbb{R})$The instance data for the class Polynomial will be the coefficients, the numbers $(a_0,…,a_N)$Provide methods that1. Evaluate the polynomial, returning $p(x)$ for any $x$2. Differentiate the polynomial, replacing the original coefficients with those of its derivative $p′$.Avoid using any `import` statements<jupyter_code>class PolyN:
    
    #use no imports
    
    def __init__(self, coeffs:list)->None:
        
        self.coeffs = coeffs
        
    def __call__(self,x:float)->float:
        
        self.expand = [(x**ix_)*cf_ for ix_,cf_ in enumerate(self.coeffs)]
        return sum(self.expand)
    
    def differentiate(self):
        
        new_coeffs = [cf_*ix_ for ix_,cf_ in enumerate(self.coeffs) if ix_>0]
        
        self.coeffs = new_coeffs
        return new_coeffs
coeffs = [1,2,3,4,5]

Polytest = PolyN(coeffs)

print(Polytest.coeffs)
print(Polytest(2))
print(Polytest.expand)
print(Polytest.differentiate())
print(Polytest.differentiate())
print(Polytest.differentiate())
print(Polytest.differentiate())
print(Polytest.differentiate())<jupyter_output>[1, 2, 3, 4, 5]
129
[1, 4, 12, 32, 80]
[2, 6, 12, 20]
[6, 24, 60]
[24, 120]
[120]
[]
{'fraction_non_alphanumeric': 0.1248366724738676, 'fraction_numerical': 0.025152439024390245, 'mean_word_length': 3.446623093681917, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 55, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '516632', 'license': 'MIT', 'n_tokens_mistral': 6584, 'n_tokens_neox': 6169, 'n_words': 2160, 'path': 'Part 3.1.ipynb', 'repo_name': 'planetnoob/QuantEcon-Python', 'size': '18368'}