Merge pull request #251 from natashawatkins/patch-7

Fix typos and prettify code in ar1_turningpts.md

Author: jstac

lectures/ar1_turningpts.md

@@ -31,9 +31,9 @@ We consider two sorts of statistics:
 
 - prospective values $y_{t+j}$ of a random process $\{y_t\}$ that is governed by the AR(1) process
 
-- sample path properties that are defined as non-linear functions of future values $\{y_{t+j}\}_{j \geq 1}$ at time $t$.
+- sample path properties that are defined as non-linear functions of future values $\{y_{t+j}\}_{j \geq 1}$ at time $t$
 
-**Sample path properties** are things like "time to next turning point" or "time to next recession"
+**Sample path properties** are things like "time to next turning point" or "time to next recession".
 
 To investigate sample path properties we'll use a simulation procedure recommended by Wecker {cite}`wecker1979predicting`.
 
@@ -43,17 +43,14 @@ Let's start with some imports.
 
 ```{code-cell} ipython3
 import numpy as np
-
 import arviz as az
 import pymc as pmc
-
 import matplotlib.pyplot as plt
 import seaborn as sns
 
 sns.set_style('white')
 colors = sns.color_palette()
 
-
 import logging
 logging.basicConfig()
 logger = logging.getLogger('pymc')
@@ -106,9 +103,9 @@ We also want to compute some  predictive distributions of "sample path statistic
 - the time until the next "recession",
 - the minimum value of $Y$ over the next 8 periods,
 - "severe recession", and
-- the time until the next turning point (positive or negative)
+- the time until the next turning point (positive or negative).
 
-To accomplish that for situations in which we are uncertain about parameter values, we shall extend the approach Wecker {cite}`wecker1979predicting` in the following way.
+To accomplish that for situations in which we are uncertain about parameter values, we shall extend Wecker's {cite}`wecker1979predicting` approach in the following way.
 
 - first simulate an initial path of length $T_0$;
 - for a given prior, draw a sample of size $N$ from the posterior joint distribution of parameters $\left(\rho,\sigma\right)$ after observing the initial path;
@@ -130,12 +127,12 @@ def AR1_simulate(rho, sigma, y0, T):
 
     # Allocate space and draw epsilons
     y = np.empty(T)
-    eps = np.random.normal(0,sigma,T)
+    eps = np.random.normal(0, sigma, T)
 
     # Initial condition and step forward
     y[0] = y0
     for t in range(1, T):
-        y[t] = rho*y[t-1] + eps[t]
+        y[t] = rho * y[t-1] + eps[t]
         
     return y
 
@@ -144,26 +141,26 @@ def plot_initial_path(initial_path):
     """
     Plot the initial path and the preceding predictive densities
     """
-    # compute .9 confidence interval]
+    # Compute .9 confidence interval]
     y0 = initial_path[-1]
-    center = np.array([rho**j*y0 for j in range(T1)])
-    vars = np.array([sigma**2*(1-rho**(2*j))/(1-rho**2) for j in range(T1)])
+    center = np.array([rho**j * y0 for j in range(T1)])
+    vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])
     y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)
     y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)
 
-    # plot
-    fig, ax = plt.subplots(1,1, figsize = (12, 6))
+    # Plot
+    fig, ax = plt.subplots(1, 1, figsize=(12, 6))
     ax.set_title("Initial Path and Predictive Densities", fontsize=15)
-    ax.plot(np.arange(-T0+1, 1), initial_path)
+    ax.plot(np.arange(-T0 + 1, 1), initial_path)
     ax.set_xlim([-T0, T1])
     ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1)
 
-    # simulate future paths
+    # Simulate future paths
     for i in range(10):
         y_future = AR1_simulate(rho, sigma, y0, T1)
         ax.plot(np.arange(T1), y_future, color='grey', alpha=.5)
     
-    # plot 90% CI
+    # Plot 90% CI
     ax.fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3, label='95% CI')
     ax.fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35, label='90% CI')
     ax.plot(np.arange(T1), center, color='red', alpha=.7, label='expected mean')
@@ -176,11 +173,11 @@ rho = 0.9
 T0, T1 = 100, 100
 y0 = 10
 
-# simulate
+# Simulate
 np.random.seed(145)
 initial_path = AR1_simulate(rho, sigma, y0, T0)
 
-# plot
+# Plot
 plot_initial_path(initial_path)
 ```
 
@@ -196,7 +193,7 @@ He called these functions "path properties" to contrast them with properties of
 
 He studied two special  prospective path properties of a given series $\{y_t\}$. 
 
-The first was **time until the next turning point**
+The first was **time until the next turning point**.
 
 * he defined a **"turning point"** to be the date of the  second of two successive declines in  $y$. 
 
@@ -227,7 +224,7 @@ W_t(\omega):= \inf \{ k\geq 1 \mid Z_{t+k}(\omega) = 1\}
 $$
 
 Wecker  {cite}`wecker1979predicting` also studied **the minimum value of $Y$ over the next 8 quarters**
-which can be defined as the random variable 
+which can be defined as the random variable.
 
 $$ 
 M_t(\omega) := \min \{ Y_{t+1}(\omega); Y_{t+2}(\omega); \dots; Y_{t+8}(\omega)\}
@@ -316,22 +313,22 @@ def draw_from_posterior(sample):
 
     with AR1_model:
         
-        #start with priors
-        rho = pmc.Uniform('rho',lower=-1.,upper=1.) #assume stable rho
+        # Start with priors
+        rho = pmc.Uniform('rho',lower=-1.,upper=1.)  # Assume stable rho
         sigma = pmc.HalfNormal('sigma', sigma = np.sqrt(10))
 
         # Expected value of y at the next period (rho * y)
-        yhat = rho*sample[:-1]
+        yhat = rho * sample[:-1]
 
         # Likelihood of the actual realization.
-        y_like = pmc.Normal('y_obs', mu = yhat, sigma=sigma, observed=sample[1:])
+        y_like = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=sample[1:])
 
     with AR1_model:
         trace = pmc.sample(10000, tune=5000)
 
     # check condition
     with AR1_model:
-        az.plot_trace(trace, figsize=(17,6))
+        az.plot_trace(trace, figsize=(17, 6))
     
     rhos = trace.posterior.rho.values.flatten()
     sigmas = trace.posterior.sigma.values.flatten()
@@ -350,7 +347,7 @@ The graphs on the left portray posterior marginal distributions.
 
 ## Calculating Sample Path Statistics
 
-Our next step is to prepare Python codeto compute our sample path statistics.
+Our next step is to prepare Python code to compute our sample path statistics.
 
 ```{code-cell} ipython3
 # define statistics
@@ -374,9 +371,9 @@ def severe_recession(omega):
     n = z.shape[0]
     
     sr = (z < -.02).astype(int)
-    indices = np.where(sr==1)[0]
+    indices = np.where(sr == 1)[0]
 
-    if len(indices)==0:
+    if len(indices) == 0:
         return T1
     else:
         return indices[0] + 1
@@ -401,8 +398,8 @@ def next_turning_point(omega):
             (omega[i+2] > omega[i+3]) and (omega[i+3] > omega[i+4])):
             T[i] = -1
     
-    up_turn = np.where(T==1)[0][0] + 1 if (1 in T) == True else T1
-    down_turn = np.where(T==-1)[0][0] + 1 if (-1 in T) == True else T1
+    up_turn = np.where(T == 1)[0][0] + 1 if (1 in T) == True else T1
+    down_turn = np.where(T == -1)[0][0] + 1 if (-1 in T) == True else T1
 
     return up_turn, down_turn
 ```
@@ -417,31 +414,31 @@ def plot_Wecker(initial_path, N, ax):
     """
     Plot the predictive distributions from "pure" Wecker's method.
     """
-    # store outcomes
+    # Store outcomes
     next_reces = np.zeros(N)
     severe_rec = np.zeros(N)
     min_vals = np.zeros(N)
     next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)
     
-    # compute .9 confidence interval]
+    # Compute .9 confidence interval]
     y0 = initial_path[-1]
-    center = np.array([rho**j*y0 for j in range(T1)])
-    vars = np.array([sigma**2*(1-rho**(2*j))/(1-rho**2) for j in range(T1)])
+    center = np.array([rho**j * y0 for j in range(T1)])
+    vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])
     y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)
     y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)
 
-    # plot
+    # Plot
     ax[0, 0].set_title("Initial path and predictive densities", fontsize=15)
-    ax[0, 0].plot(np.arange(-T0+1, 1), initial_path)
+    ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)
     ax[0, 0].set_xlim([-T0, T1])
     ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)
 
-    # plot 90% CI
+    # Plot 90% CI
     ax[0, 0].fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3)
     ax[0, 0].fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35)
     ax[0, 0].plot(np.arange(T1), center, color='red', alpha=.7)
 
-    # simulate future paths
+    # Simulate future paths
     for n in range(N):
         sim_path = AR1_simulate(rho, sigma, initial_path[-1], T1)
         next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))
@@ -452,7 +449,7 @@ def plot_Wecker(initial_path, N, ax):
         if n%(N/10) == 0:
             ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)
             
-    # return next_up_turn, next_down_turn
+    # Return next_up_turn, next_down_turn
     sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.8, label='True parameters')
     ax[0, 1].set_title("Predictive distribution of time until the next recession", fontsize=13)
 
@@ -478,42 +475,42 @@ plt.show()
 Now we apply we apply our  "extended" Wecker method based on  predictive densities of $y$ defined by
 {eq}`ar1-tp-eq4` that acknowledge posterior uncertainty in the parameters $\rho, \sigma$.
 
-To approximate  the intergration on the right side of {eq}`ar1-tp-eq4`, we  repeately draw parameters from the joint posterior distribution each time we simulate a sequence of future values from model {eq}`ar1-tp-eq1`.
+To approximate  the intergration on the right side of {eq}`ar1-tp-eq4`, we  repeatedly draw parameters from the joint posterior distribution each time we simulate a sequence of future values from model {eq}`ar1-tp-eq1`.
 
 ```{code-cell} ipython3
 def plot_extended_Wecker(post_samples, initial_path, N, ax):
     """
     Plot the extended Wecker's predictive distribution
     """
-    # select a sample
-    index = np.random.choice(np.arange(len(post_samples['rho'])), N+1, replace=False)
+    # Select a sample
+    index = np.random.choice(np.arange(len(post_samples['rho'])), N + 1, replace=False)
     rho_sample = post_samples['rho'][index]
     sigma_sample = post_samples['sigma'][index]
 
-    # store outcomes
+    # Store outcomes
     next_reces = np.zeros(N)
     severe_rec = np.zeros(N)
     min_vals = np.zeros(N)
     next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)
 
-    # plot
+    # Plot
     ax[0, 0].set_title("Initial path and future paths simulated from posterior draws", fontsize=15)
-    ax[0, 0].plot(np.arange(-T0+1, 1), initial_path)
+    ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)
     ax[0, 0].set_xlim([-T0, T1])
     ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)
 
-    # simulate future paths
+    # Simulate future paths
     for n in range(N):
         sim_path = AR1_simulate(rho_sample[n], sigma_sample[n], initial_path[-1], T1)
         next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))
         severe_rec[n] = severe_recession(sim_path)
         min_vals[n] = minimum_value(sim_path)
         next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path)
 
-        if n%(N/10) == 0:
+        if n % (N / 10) == 0:
             ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)
         
-    # return next_up_turn, next_down_turn
+    # Return next_up_turn, next_down_turn
     sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.6, color=colors[1], label='Sampling from posterior')
     ax[0, 1].set_title("Predictive distribution of time until the next recession", fontsize=13)
 
@@ -529,19 +526,19 @@ def plot_extended_Wecker(post_samples, initial_path, N, ax):
     sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.6, color=colors[1], label='Sampling from posterior')
     ax[2, 1].set_title("Predictive distribution of time until the next negative turn", fontsize=13)
 
-fig, ax = plt.subplots(3, 2, figsize=(15,12))
+fig, ax = plt.subplots(3, 2, figsize=(15, 12))
 plot_extended_Wecker(post_samples, initial_path, 1000, ax)
 plt.show()
 ```
 
 ## Comparison
 
-Finally, we plot both the original Wecker method and the extended method with parameter values drawn from the posterior together to compare the differnces that emerge from pretending to know parameter values when they are actually uncertain.  
+Finally, we plot both the original Wecker method and the extended method with parameter values drawn from the posterior together to compare the differences that emerge from pretending to know parameter values when they are actually uncertain.  
 
 ```{code-cell} ipython3
 fig, ax = plt.subplots(3, 2, figsize=(15,12))
 plot_Wecker(initial_path, 1000, ax)
-ax[0,0].clear()
+ax[0, 0].clear()
 plot_extended_Wecker(post_samples, initial_path, 1000, ax)
 plt.legend()
 plt.show()