From fc930ae85b94755daff6b18a043ae03bd4175559 Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Tue, 26 Aug 2025 17:24:47 +0800 Subject: [PATCH 01/16] Update ar1_turningpts.md --- lectures/ar1_turningpts.md | 151 +++++++++++++++++++++---------------- 1 file changed, 85 insertions(+), 66 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 3aa55a9df..f66d8ade5 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -48,6 +48,15 @@ import pymc as pmc import matplotlib.pyplot as plt import seaborn as sns +# jax +import jax.random as random +import jax.numpy as jnp + +# numpyro +import numpyro +import numpyro.distributions as dist +from numpyro.infer import MCMC, NUTS + sns.set_style('white') colors = sns.color_palette() @@ -91,12 +100,16 @@ We form this predictive distribution by integrating {eq}`ar1-tp-eq3` with respec that conditions on an observed history $y^t = \{y_s\}_{s=0}^t$: $$ -f(y_{t+j} | y^t) = \int f(y_{t+j} | y_{t}; \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma +\begin{aligned} +f(y_{t+j} | y^t) += \int f(y_{t+j} | y^t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma += \int f(y_{t+j} | y_t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma +\end{aligned} $$ (ar1-tp-eq4) -Predictive distribution {eq}`ar1-tp-eq3` assumes that parameters $(\rho,\sigma)$ are known. +Predictive distribution {eq}`ar1-tp-eq3` assumes that parameters $(\rho,\sigma)$ are known. -Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. +Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows that $\{y_t\}$ is a AR(1) process when $\(\rho, \sigma\)$ are given. We also want to compute some predictive distributions of "sample path statistics" that might include, for example @@ -123,8 +136,22 @@ we'll plot $.9$ and $.95$ coverage intervals using conditional distribution We'll also plot a bunch of samples of sequences of future values and watch where they fall relative to the coverage interval. ```{code-cell} ipython3 -def AR1_simulate(rho, sigma, y0, T): +class AR1(NamedTuple): + """Define a AR1 process + + Parameters + --------- + rho: float + coefficient + sigma: float + standard error of the error term + """ + rho:float + sigma:float + +def AR1_simulate(rho, sigma, y0, T): + """Generate a realization of the AR1 process given parameters and length""" # Allocate space and draw epsilons y = np.empty(T) eps = np.random.normal(0, sigma, T) @@ -134,42 +161,63 @@ def AR1_simulate(rho, sigma, y0, T): for t in range(1, T): y[t] = rho * y[t-1] + eps[t] - return y + return jnp.asarray(y) -def plot_initial_path(initial_path): - """ - Plot the initial path and the preceding predictive densities - """ - # 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)]) - 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) +def plot_initial_path(ar1, initial_path, predict_length): + """Plot the initial path and the preceding predictive densities""" + rho, sigma = ar1.rho, ar1.sigma + T0 = len(initial_path) + T1 = predict_length + + # Compute moments and confidence intervals + y_T0 = initial_path[-1] + center = jnp.array([rho**j * y_T0 for j in range(T1)]) + vars = jnp.array( + [sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)] + ) + + y_upper_c95 = center + 1.96 * np.sqrt(vars) + y_lower_c95 = center - 1.96 * np.sqrt(vars) + + y_upper_c90 = center + 1.65 * np.sqrt(vars) + y_lower_c90 = center - 1.65 * np.sqrt(vars) # Plot - fig, ax = plt.subplots(1, 1, figsize=(12, 6)) - ax.set_title("Initial Path and Predictive Densities", fontsize=15) + fig, ax = plt.subplots(1, 1) + #ax.set_title("Initial Path and Predictive Densities", fontsize=15) 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 10 future paths for i in range(10): - y_future = AR1_simulate(rho, sigma, y0, T1) + y_future = AR1_simulate(rho, sigma, y_T0, T1) ax.plot(np.arange(T1), y_future, color='grey', alpha=.5) # 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') + ax.fill_between( + np.arange(T1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' + ) + ax.fill_between( + np.arange(T1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' + ) + ax.plot(np.arange(T1), center, color='red', alpha=.7, label='expectation') ax.legend(fontsize=12) plt.show() +``` - -sigma = 1 +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Initial path and predictive densities + name: fig_path +--- +ar1 = AR1(rho=0.9, sigma=1) rho = 0.9 +simga = 1 T0, T1 = 100, 100 y0 = 10 @@ -178,7 +226,7 @@ np.random.seed(145) initial_path = AR1_simulate(rho, sigma, y0, T0) # Plot -plot_initial_path(initial_path) +plot_initial_path(ar1, initial_path, T1) ``` As functions of forecast horizon, the coverage intervals have shapes like those described in @@ -199,25 +247,17 @@ The first was **time until the next turning point**. To examine this statistic, let $Z$ be an indicator process - - - - $$ -Z_t(Y(\omega)) := +Z_t(\omega) := \begin{cases} \ 1 & \text{if } Y_t(\omega)< Y_{t-1}(\omega)< Y_{t-2}(\omega) \geq Y_{t-3}(\omega) \\ -0 & \text{otherwise} +\ 0 & \text{otherwise} \end{cases} $$ -Then the random variable **time until the next turning point** is defined as the following **stopping time** with respect to $Z$: +Here $\omega \in Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow R$ gives $y_t$ according to $\omega$ and the AR(1) process. + +Then the random variable **time until the next turning point** is defined as the following *stopping time* with respect to $Z$: $$ W_t(\omega):= \inf \{ k\geq 1 \mid Z_{t+k}(\omega) = 1\} @@ -233,50 +273,29 @@ $$ It is interesting to study yet another possible concept of a **turning point**. Thus, let - - $$ -T_t(Y(\omega)) := +T_t(\omega) := \begin{cases} \ 1 & \text{if } Y_{t-2}(\omega)> Y_{t-1}(\omega) > Y_{t}(\omega) \ \text{and } \ Y_{t}(\omega) < Y_{t+1}(\omega) < Y_{t+2}(\omega) \\ \ -1 & \text{if } Y_{t-2}(\omega)< Y_{t-1}(\omega) < Y_{t}(\omega) \ \text{and } \ Y_{t}(\omega) > Y_{t+1}(\omega) > Y_{t+2}(\omega) \\ -0 & \text{otherwise} +\ 0 & \text{otherwise} \end{cases} $$ - - Define a **positive turning point today or tomorrow** statistic as - - - $$ P_t(\omega) := \begin{cases} \ 1 & \text{if } T_t(\omega)=1 \ \text{or} \ T_{t+1}(\omega)=1 \\ -0 & \text{otherwise} +\ 0 & \text{otherwise} \end{cases} $$ This is designed to express the event -- ``after one or two decrease(s), $Y$ will grow for two consecutive quarters'' +- "after one or two decrease(s), $Y$ will grow for two consecutive quarters" Following {cite}`wecker1979predicting`, we can use simulations to calculate probabilities of $P_t$ and $N_t$ for each period $t$. @@ -286,15 +305,15 @@ The procedure consists of the following steps: * index a sample path by $\omega_i$ -* for a given date $t$, simulate $I$ sample paths of length $N$ +* from a given date $t$, simulate $I$ sample paths of length $N$ $$ Y(\omega_i) = \left\{ Y_{t+1}(\omega_i), Y_{t+2}(\omega_i), \dots, Y_{t+N}(\omega_i)\right\}_{i=1}^I $$ -* for each path $\omega_i$, compute the associated value of $W_t(\omega_i), W_{t+1}(\omega_i), \dots$ +* for each path $\omega_i$, compute the associated value of $W_t(\omega_i), W_{t+1}(\omega_i), \dots , W_{t+N}$ -* consider the sets $\{W_t(\omega_i)\}^{T}_{i=1}, \ \{W_{t+1}(\omega_i)\}^{T}_{i=1}, \ \dots, \ \{W_{t+N}(\omega_i)\}^{T}_{i=1}$ as samples from the predictive distributions $f(W_{t+1} \mid \mathcal y_t, \dots)$, $f(W_{t+2} \mid y_t, y_{t-1}, \dots)$, $\dots$, $f(W_{t+N} \mid y_t, y_{t-1}, \dots)$. +* consider the sets $\{W_t(\omega_i)\}^{I}_{i=1}, \ \{W_{t+1}(\omega_i)\}^{I}_{i=1}, \ \dots, \ \{W_{t+N}(\omega_i)\}^{I}_{i=1}$ as samples from the predictive distributions $f(W_{t+1} \mid y_t, y_{t-1}, \dots , y_0)$, $f(W_{t+2} \mid y_t, y_{t-1}, \dots , y_0)$, $\dots$, $f(W_{t+N} \mid y_t, y_{t-1}, \dots , y_0)$. ## Using Simulations to Approximate a Posterior Distribution From 9700a3aa34f72a643431b53f1c715793e9617d0b Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Wed, 27 Aug 2025 17:30:57 +0800 Subject: [PATCH 02/16] Update ar1_turningpts.md --- lectures/ar1_turningpts.md | 179 +++++++++++++++++++------------------ 1 file changed, 94 insertions(+), 85 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index f66d8ade5..4b121020f 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -59,11 +59,7 @@ from numpyro.infer import MCMC, NUTS sns.set_style('white') colors = sns.color_palette() - -import logging -logging.basicConfig() -logger = logging.getLogger('pymc') -logger.setLevel(logging.CRITICAL) +key = random.PRNGKey(0) ``` ## A Univariate First-Order Autoregressive Process @@ -150,21 +146,21 @@ class AR1(NamedTuple): sigma:float -def AR1_simulate(rho, sigma, y0, T): +def AR1_simulate(ar1: AR1, y0, T, subkey): """Generate a realization of the AR1 process given parameters and length""" # Allocate space and draw epsilons - y = np.empty(T) - eps = np.random.normal(0, sigma, T) + y = jnp.empty(T) + eps = ar1.sigma * random.normal(subkey, (T,)) # Initial condition and step forward - y[0] = y0 + y.at[0].set(y0) for t in range(1, T): - y[t] = rho * y[t-1] + eps[t] + y.at[t].set(ar1.rho * y[t-1] + eps[t]) return jnp.asarray(y) -def plot_initial_path(ar1, initial_path, predict_length): +def plot_path(ar1, initial_path, predict_length): """Plot the initial path and the preceding predictive densities""" rho, sigma = ar1.rho, ar1.sigma T0 = len(initial_path) @@ -177,32 +173,33 @@ def plot_initial_path(ar1, initial_path, predict_length): [sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)] ) - y_upper_c95 = center + 1.96 * np.sqrt(vars) - y_lower_c95 = center - 1.96 * np.sqrt(vars) + y_upper_c95 = center + 1.96 * jnp.sqrt(vars) + y_lower_c95 = center - 1.96 * jnp.sqrt(vars) - y_upper_c90 = center + 1.65 * np.sqrt(vars) - y_lower_c90 = center - 1.65 * np.sqrt(vars) + y_upper_c90 = center + 1.65 * jnp.sqrt(vars) + y_lower_c90 = center - 1.65 * jnp.sqrt(vars) # Plot fig, ax = plt.subplots(1, 1) #ax.set_title("Initial Path and Predictive Densities", fontsize=15) - ax.plot(np.arange(-T0 + 1, 1), initial_path) + ax.plot(jnp.arange(-T0 + 1, 1), initial_path) ax.set_xlim([-T0, T1]) ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1) # Simulate 10 future paths + subkeys = random.split(key, num=10) for i in range(10): - y_future = AR1_simulate(rho, sigma, y_T0, T1) - ax.plot(np.arange(T1), y_future, color='grey', alpha=.5) + y_future = AR1_simulate(ar1, y_T0, T1, subkeys[i]) + ax.plot(jnp.arange(T1), y_future, color='grey', alpha=.5) # Plot 90% CI ax.fill_between( - np.arange(T1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' + jnp.arange(T1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' ) ax.fill_between( - np.arange(T1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' + jnp.arange(T1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' ) - ax.plot(np.arange(T1), center, color='red', alpha=.7, label='expectation') + ax.plot(jnp.arange(T1), center, color='red', alpha=.7, label='expectation') ax.legend(fontsize=12) plt.show() ``` @@ -216,17 +213,14 @@ mystnb: name: fig_path --- ar1 = AR1(rho=0.9, sigma=1) -rho = 0.9 -simga = 1 T0, T1 = 100, 100 y0 = 10 # Simulate -np.random.seed(145) -initial_path = AR1_simulate(rho, sigma, y0, T0) +initial_path = AR1_simulate(ar1, y0, T0) # Plot -plot_initial_path(ar1, initial_path, T1) +plot_path(ar1, initial_path, T1) ``` As functions of forecast horizon, the coverage intervals have shapes like those described in @@ -323,110 +317,122 @@ The next code cells use `pymc` to compute the time $t$ posterior distribution of Note that in defining the likelihood function, we choose to condition on the initial value $y_0$. ```{code-cell} ipython3 -def draw_from_posterior(sample): - """ - Draw a sample of size N from the posterior distribution. - """ - - AR1_model = pmc.Model() - - with AR1_model: - - # 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] +--- +mystnb: + figure: + caption: | + Posterior distributions (rho, sigma) + name: fig_post +--- +def draw_from_posterior(sample, size=1e5, bins=20, dis_plot=1): + """Draw a sample of size from the posterior distribution.""" + # Start with priors + rho = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable rho + sigma = numpyro.sample('sigma', dis.HalfNormal(sqrt(10))) + + # Define likelihood recursively + for t in range(1, len(sample)) + # Expectation of y_t + mu = rho * sample[t-1] # Likelihood of the actual realization. - y_like = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=sample[1:]) - - with AR1_model: - trace = pmc.sample(10000, tune=5000) + numpyro.sample(f'y_{t}', dist.Normal(mu, sigma), obs=sample[t]) - # check condition - with AR1_model: - az.plot_trace(trace, figsize=(17, 6)) - - rhos = trace.posterior.rho.values.flatten() - sigmas = trace.posterior.sigma.values.flatten() + # Compute posterior distribution of parameters + nuts_kernel = NUTS(model) + mcmc = MCMC( + nuts_kernek, + num_warmup=5000, + num_sample=size, + progress_bar=False) + mcmc.run(random.PRNGKey(0), sample=sample) post_sample = { - 'rho': rhos, - 'sigma': sigmas + 'rho': mcmc.get_samples()['rho'], + 'sigma': mcmc.get_samples()['sigma'] } - + + if dis_plot==1: + sns.displot( + post_sample, kde=True, stat="density", bins=bins, height=5, aspect=1.5 + ) + return post_sample post_samples = draw_from_posterior(initial_path) ``` -The graphs on the left portray posterior marginal distributions. +The graphs above portray posterior distributions. ## Calculating Sample Path Statistics Our next step is to prepare Python code to compute our sample path statistics. +These statistics were originally defined as random variables with respect to $\omega$, but here we use $\{Y_t\}$ as the argument because $\omega$ is implicit. + +Also, these two kinds of definitions are equivalent because $\omega$ determins path statistics only through $\{Y_t\}$ + ```{code-cell} ipython3 -# define statistics -def next_recession(omega): - n = omega.shape[0] - 3 - z = np.zeros(n, dtype=int) +# Define statistics as functions of y, the realized path values. +def next_recession(y): + n = y.shape[0] - 3 + z = jnp.zeros(n, dtype=int) for i in range(n): - z[i] = int(omega[i] <= omega[i+1] and omega[i+1] > omega[i+2] and omega[i+2] > omega[i+3]) + z.at[i].set(int(y[i] <= y[i+1] and y[i+1] > y[i+2] and y[i+2] > y[i+3])) - if np.any(z) == False: + if jnp.any(z) == False: return 500 else: - return np.where(z==1)[0][0] + 1 + return jnp.where(z == 1)[0][0] + 1 -def minimum_value(omega): - return min(omega[:8]) -def severe_recession(omega): - z = np.diff(omega) +def next_severe_recession(y): + z = jnp.diff(y) n = z.shape[0] sr = (z < -.02).astype(int) - indices = np.where(sr == 1)[0] + indices = jnp.where(sr == 1)[0] if len(indices) == 0: return T1 else: return indices[0] + 1 -def next_turning_point(omega): + +def minimum_value_8q(y): + return jnp.min(y[:8]) + + +def next_turning_point(y): """ - Suppose that omega is of length 6 + Suppose that y is of length 6 y_{t-2}, y_{t-1}, y_{t}, y_{t+1}, y_{t+2}, y_{t+3} that is sufficient for determining the value of P/N """ - n = np.asarray(omega).shape[0] - 4 - T = np.zeros(n, dtype=int) + n = jnp.asarray(y).shape[0] - 4 + T = jnp.zeros(n, dtype=int) for i in range(n): - if ((omega[i] > omega[i+1]) and (omega[i+1] > omega[i+2]) and - (omega[i+2] < omega[i+3]) and (omega[i+3] < omega[i+4])): + if ((y[i] > y[i+1]) and (y[i+1] > y[i+2]) and + (y[i+2] < y[i+3]) and (y[i+3] < y[i+4])): T[i] = 1 - elif ((omega[i] < omega[i+1]) and (omega[i+1] < omega[i+2]) and - (omega[i+2] > omega[i+3]) and (omega[i+3] > omega[i+4])): + elif ((y[i] < y[i+1]) and (y[i+1] < y[i+2]) and + (y[i+2] > y[i+3]) and (y[i+3] > y[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 = jnp.where(T == 1)[0][0] + 1 if (1 in T) == True else T1 + down_turn = jnp.where(T == -1)[0][0] + 1 if (-1 in T) == True else T1 return up_turn, down_turn ``` ## Original Wecker Method -Now we apply Wecker's original method by simulating future paths and compute predictive distributions, conditioning -on the true parameters associated with the data-generating model. +Now we apply Wecker's original method by simulating future paths and compute predictive distributions, conditioning on the true parameters associated with the data-generating model. ```{code-cell} ipython3 def plot_Wecker(initial_path, N, ax): @@ -458,11 +464,12 @@ def plot_Wecker(initial_path, N, ax): ax[0, 0].plot(np.arange(T1), center, color='red', alpha=.7) # Simulate future paths + subkeys = random.split(key, num=N) for n in range(N): - sim_path = AR1_simulate(rho, sigma, initial_path[-1], T1) + sim_path = AR1_simulate(ar1, initial_path[-1], T1, subkeys[n]) 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) + severe_rec[n] = next_severe_recession(sim_path) + min_vals[n] = minimum_value_8q(sim_path) next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) if n%(N/10) == 0: @@ -519,11 +526,13 @@ def plot_extended_Wecker(post_samples, initial_path, N, ax): ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1) # Simulate future paths + subkeys = random.split(key, num=N) for n in range(N): - sim_path = AR1_simulate(rho_sample[n], sigma_sample[n], initial_path[-1], T1) + ar1_n = AR1(rho=rho_sample[n], sigma=sigma_sample[n]) + sim_path = AR1_simulate(ar1_n, initial_path[-1], T1, subkeys[n]) 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) + severe_rec[n] = next_severe_recession(sim_path) + min_vals[n] = minimum_value_8q(sim_path) next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) if n % (N / 10) == 0: From 2150930b72d463301b1342b31030afc4202e8822 Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Fri, 29 Aug 2025 16:45:32 +0800 Subject: [PATCH 03/16] Update ar1_turningpts.md --- lectures/ar1_turningpts.md | 297 +++++++++++++++++++++++++++---------- 1 file changed, 216 insertions(+), 81 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 4b121020f..98847ae84 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -133,21 +133,49 @@ We'll also plot a bunch of samples of sequences of future values and watch where ```{code-cell} ipython3 class AR1(NamedTuple): - """Define a AR1 process + """ + Represents a univariate first-order autoregressive (AR(1)) process. Parameters - --------- - rho: float - coefficient - sigma: float - standard error of the error term + ---------- + rho : float + Autoregressive coefficient, must satisfy |rho| < 1 for stationarity. + sigma : float + Standard deviation of the error term. + y0 : float + Initial value of the process at time t=0. + T0 : int, optional + Length of the initial observed path (default is 100). + T1 : int, optional + Length of the future path to simulate (default is 100). """ rho:float sigma:float + y0: float + T0: int=100 + T1: int=100 + +def AR1_simulate(ar1: AR1, T, subkey): + """ + Simulate a realization of the AR(1) process for a given number of periods. -def AR1_simulate(ar1: AR1, y0, T, subkey): - """Generate a realization of the AR1 process given parameters and length""" + Parameters + ---------- + ar1 : AR1 + AR1 named tuple containing process parameters (rho, sigma, y0). + T : int + Number of time steps to simulate. + subkey : jax.random.PRNGKey + JAX random key for generating random noise. + + Returns + ------- + y : jax.numpy.ndarray + Simulated path of the AR(1) process of length T. + """ + # ...function code... + rho, sigma, y0 = ar.rho, ar1.sigma, ar1.y0 # Allocate space and draw epsilons y = jnp.empty(T) eps = ar1.sigma * random.normal(subkey, (T,)) @@ -160,11 +188,9 @@ def AR1_simulate(ar1: AR1, y0, T, subkey): return jnp.asarray(y) -def plot_path(ar1, initial_path, predict_length): +def plot_path(ar1, initial_path, key, ax): """Plot the initial path and the preceding predictive densities""" - rho, sigma = ar1.rho, ar1.sigma - T0 = len(initial_path) - T1 = predict_length + rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 # Compute moments and confidence intervals y_T0 = initial_path[-1] @@ -180,8 +206,6 @@ def plot_path(ar1, initial_path, predict_length): y_lower_c90 = center - 1.65 * jnp.sqrt(vars) # Plot - fig, ax = plt.subplots(1, 1) - #ax.set_title("Initial Path and Predictive Densities", fontsize=15) ax.plot(jnp.arange(-T0 + 1, 1), initial_path) ax.set_xlim([-T0, T1]) ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1) @@ -189,7 +213,8 @@ def plot_path(ar1, initial_path, predict_length): # Simulate 10 future paths subkeys = random.split(key, num=10) for i in range(10): - y_future = AR1_simulate(ar1, y_T0, T1, subkeys[i]) + ar1_n = AR1(rho=rho, sigma=sigma, y0=y_T0) + y_future = AR1_simulate(ar1, T1, subkeys[i]) ax.plot(jnp.arange(T1), y_future, color='grey', alpha=.5) # Plot 90% CI @@ -201,7 +226,6 @@ def plot_path(ar1, initial_path, predict_length): ) ax.plot(jnp.arange(T1), center, color='red', alpha=.7, label='expectation') ax.legend(fontsize=12) - plt.show() ``` ```{code-cell} ipython3 @@ -212,15 +236,15 @@ mystnb: Initial path and predictive densities name: fig_path --- -ar1 = AR1(rho=0.9, sigma=1) -T0, T1 = 100, 100 -y0 = 10 +ar1 = AR1(rho=0.9, sigma=1, y0=10, T0=100, T1=100) # Simulate -initial_path = AR1_simulate(ar1, y0, T0) +initial_path = AR1_simulate(ar1, ar1.T0, key) # Plot -plot_path(ar1, initial_path, T1) +fig, ax = plt.subplots(1, 1) +plot_path(ar1, initial_path, key, ax) +plt.show() ``` As functions of forecast horizon, the coverage intervals have shapes like those described in @@ -356,9 +380,10 @@ def draw_from_posterior(sample, size=1e5, bins=20, dis_plot=1): sns.displot( post_sample, kde=True, stat="density", bins=bins, height=5, aspect=1.5 ) - + # TODO:Check plot and return return post_sample + post_samples = draw_from_posterior(initial_path) ``` @@ -435,64 +460,118 @@ def next_turning_point(y): Now we apply Wecker's original method by simulating future paths and compute predictive distributions, conditioning on the true parameters associated with the data-generating model. ```{code-cell} ipython3 -def plot_Wecker(initial_path, N, ax): +--- +mystnb: + figure: + caption: | + Distributions of statistics by Wecker's method + name: fig_wecker +--- +def plot_Wecker(ar1: AR1, initial_path, N, ax): """ Plot the predictive distributions from "pure" Wecker's method. + + Parameters + ---------- + ar1 : AR1 + An AR1 named tuple containing the process parameters (rho, sigma, T0, T1). + initial_path : array-like + The initial observed path of the AR(1) process. + N : int + Number of future sample paths to simulate for predictive distributions. """ + rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 # 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) + next_reces = jnp.zeros(N) + severe_rec = jnp.zeros(N) + min_vals = jnp.zeros(N) + next_up_turn, next_down_turn = jnp.zeros(N), jnp.zeros(N) - # 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)]) - 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) + # Compute confidence intervals + y_T0 = initial_path[-1] + center = jnp.array([rho**j * yT0 for j in range(T1)]) + vars = jnp.array( + [sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)] + ) + y_upper_c95 = center + 1.96 * jnp.sqrt(vars) + y_lower_c95 = center - 1.96 * jnp.sqrt(vars) + + y_upper_c90 = center + 1.65 * jnp.sqrt(vars) + y_lower_c90 = center - 1.65 * jnp.sqrt(vars) # 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(jnp.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 - 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) + ax[0, 0].fill_between(jnp.arange(T1), y_upper_c95, y_lower_c95, alpha=.3) + ax[0, 0].fill_between(jnp.arange(T1), y_upper_c90, y_lower_c90, alpha=.35) + ax[0, 0].plot(jnp.arange(T1), center, color='red', alpha=.7) # Simulate future paths subkeys = random.split(key, num=N) for n in range(N): - sim_path = AR1_simulate(ar1, initial_path[-1], T1, subkeys[n]) - next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path])) + ar1_n = AR1(rho=rho, sigma=sigma, y0=y_T0) + sim_path = AR1_simulate(ar1_n, T1, subkeys[n]) + next_reces[n] = next_recession(jnp.hstack([initial_path[-3:-1], sim_path])) severe_rec[n] = next_severe_recession(sim_path) min_vals[n] = minimum_value_8q(sim_path) next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) if n%(N/10) == 0: - ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1) + ax[0, 0].plot(jnp.arange(T1), sim_path, color='gray', alpha=.3, lw=1) # 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) + sns.histplot( + next_reces, kde=True, stat='density', + ax=ax[0, 1], alpha=.8, label='True parameters' + ) + ax[0, 1].set_title( + "Predictive distribution (time until the next recession)", + fontsize=13 + ) + + sns.histplot( + severe_rec, kde=False, stat='density', + ax=ax[1, 0], binwidth=0.9, alpha=.7, label='True parameters' + ) + ax[1, 0].set_title( + r"Predictive distribution (stopping time of growth$<-2\%$)", + fontsize=13 + ) - sns.histplot(severe_rec, kde=False, stat='density', ax=ax[1, 0], binwidth=0.9, alpha=.7, label='True parameters') - ax[1, 0].set_title(r"Predictive distribution of stopping time of growth$<-2\%$", fontsize=13) + sns.histplot( + min_vals, kde=True, stat='density', + ax=ax[1, 1], alpha=.8, label='True parameters' + ) + ax[1, 1].set_title( + "Predictive distribution (minimum value in next 8 periods)", + fontsize=13 + ) - sns.histplot(min_vals, kde=True, stat='density', ax=ax[1, 1], alpha=.8, label='True parameters') - ax[1, 1].set_title("Predictive distribution of minimum value in the next 8 periods", fontsize=13) + sns.histplot( + next_up_turn, kde=True, stat='density', + ax=ax[2, 0], alpha=.8, label='True parameters' + ) + ax[2, 0].set_title( + "Predictive distribution (time until the next positive turn)", + fontsize=13 + ) - sns.histplot(next_up_turn, kde=True, stat='density', ax=ax[2, 0], alpha=.8, label='True parameters') - ax[2, 0].set_title("Predictive distribution of time until the next positive turn", fontsize=13) + sns.histplot( + next_down_turn, kde=True, stat='density', + ax=ax[2, 1], alpha=.8, label='True parameters' + ) + ax[2, 1].set_title( + "Predictive distribution (time until the next negative turn)", + fontsize=13 + ) - sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.8, label='True parameters') - 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)) -plot_Wecker(initial_path, 1000, ax) +fig, ax = plt.subplots(3, 2) +plot_Wecker(ar1, initial_path, 1000) plt.show() ``` @@ -504,58 +583,107 @@ Now we apply we apply our "extended" Wecker method based on predictive densiti 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): +--- +mystnb: + figure: + caption: | + Distributions of statistics by extended Wecker's method + name: fig_extend_wecker +--- +def plot_extended_Wecker(ar1, post_samples, initial_path, N, ax): """ Plot the extended Wecker's predictive distribution """ + rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 # Select a sample - index = np.random.choice(np.arange(len(post_samples['rho'])), N + 1, replace=False) + index = random.choice( + key, jnp.arange(len(post_samples['rho'])), (N + 1,), replace=False + ) rho_sample = post_samples['rho'][index] sigma_sample = post_samples['sigma'][index] # 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) + next_reces = jnp.zeros(N) + severe_rec = jnp.zeros(N) + min_vals = jnp.zeros(N) + next_up_turn, next_down_turn = jnp.zeros(N), jnp.zeros(N) # 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].set_title( + "Initial path and future paths simulated from posterior draws", + fontsize=15 + ) + ax[0, 0].plot(jnp.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 subkeys = random.split(key, num=N) for n in range(N): - ar1_n = AR1(rho=rho_sample[n], sigma=sigma_sample[n]) - sim_path = AR1_simulate(ar1_n, initial_path[-1], T1, subkeys[n]) - next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path])) + ar1_n = AR1( + rho=rho_sample[n], sigma=sigma_sample[n], y0= initial_path[-1] + ) + sim_path = AR1_simulate(ar1_n, T1, subkeys[n]) + next_reces[n] = next_recession( + jnp.hstack([initial_path[-3:-1], sim_path]) + ) severe_rec[n] = next_severe_recession(sim_path) min_vals[n] = minimum_value_8q(sim_path) next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) if n % (N / 10) == 0: - ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1) + ax[0, 0].plot(jnp.arange(T1), sim_path, color='gray', alpha=.3, lw=1) # 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) - - sns.histplot(severe_rec, kde=False, stat='density', ax=ax[1, 0], binwidth=.9, alpha=.6, color=colors[1], label='Sampling from posterior') - ax[1, 0].set_title(r"Predictive distribution of stopping time of growth$<-2\%$", fontsize=13) - - sns.histplot(min_vals, kde=True, stat='density', ax=ax[1, 1], alpha=.6, color=colors[1], label='Sampling from posterior') - ax[1, 1].set_title("Predictive distribution of minimum value in the next 8 periods", fontsize=13) + 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 (time until the next recession)", + fontsize=13) + + sns.histplot( + severe_rec, kde=False, stat='density', + ax=ax[1, 0], binwidth=.9, alpha=.6, color=colors[1], + label='Sampling from posterior' + ) + ax[1, 0].set_title( + r"Predictive distribution (stopping time of growth$<-2\%$)", + fontsize=13 + ) - sns.histplot(next_up_turn, kde=True, stat='density', ax=ax[2, 0], alpha=.6, color=colors[1], label='Sampling from posterior') - ax[2, 0].set_title("Predictive distribution of time until the next positive turn", fontsize=13) + sns.histplot( + min_vals, kde=True, stat='density', + ax=ax[1, 1], alpha=.6, color=colors[1], + label='Sampling from posterior' + ) + ax[1, 1].set_title( + "Predictive distribution (minimum value in the next 8 periods)", + fontsize=13) + + sns.histplot( + next_up_turn, kde=True, stat='density', + ax=ax[2, 0], alpha=.6, color=colors[1], + label='Sampling from posterior') + ax[2, 0].set_title( + "Predictive distribution of time until the next positive turn", + fontsize=13 + ) - 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) + 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)) -plot_extended_Wecker(post_samples, initial_path, 1000, ax) +fig, ax = plt.subplots(3, 2) +plot_extended_Wecker(ar1, post_samples, initial_path, 1000, ax) plt.show() ``` @@ -564,10 +692,17 @@ plt.show() 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) +--- +mystnb: + figure: + caption: | + Comparison between two methods + name: fig_wecker +--- +fig, ax = plt.subplots(3, 2) +plot_Wecker(ar1, initial_path, 1000, ax) ax[0, 0].clear() -plot_extended_Wecker(post_samples, initial_path, 1000, ax) +plot_extended_Wecker(ar1, post_samples, initial_path, 1000, ax) plt.legend() plt.show() ``` From 9d351da60acb33ce31b5682330af0eb767a03764 Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Sun, 31 Aug 2025 16:18:25 +0800 Subject: [PATCH 04/16] Update ar1_turningpts.md --- lectures/ar1_turningpts.md | 563 ++++++++++++++++++------------------- 1 file changed, 266 insertions(+), 297 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 98847ae84..2dac59d9a 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -1,10 +1,10 @@ --- jupytext: text_representation: - extension: .myst + extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.13.8 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -13,10 +13,13 @@ kernelspec: # Forecasting an AR(1) Process +```{include} _admonition/gpu.md +``` + ```{code-cell} ipython3 :tags: [hide-output] -!pip install arviz pymc +!pip install numpyro jax ``` This lecture describes methods for forecasting statistics that are functions of future values of a univariate autogressive process. @@ -42,15 +45,15 @@ To acknowledge uncertainty about parameters, we'll deploy `pymc` to construct a 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 +from typing import NamedTuple # jax +import jax import jax.random as random import jax.numpy as jnp +from jax import lax # numpyro import numpyro @@ -98,14 +101,14 @@ that conditions on an observed history $y^t = \{y_s\}_{s=0}^t$: $$ \begin{aligned} f(y_{t+j} | y^t) -= \int f(y_{t+j} | y^t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma -= \int f(y_{t+j} | y_t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma +&= \int f(y_{t+j} | y^t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma\\ +&= \int f(y_{t+j} | y_t, \rho, \sigma) \pi_t(\rho,\sigma | y^t ) d \rho d \sigma \end{aligned} $$ (ar1-tp-eq4) Predictive distribution {eq}`ar1-tp-eq3` assumes that parameters $(\rho,\sigma)$ are known. -Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows that $\{y_t\}$ is a AR(1) process when $\(\rho, \sigma\)$ are given. +Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows that $\{y_t\}$ is a AR(1) process when $(\rho, \sigma)$ are given. We also want to compute some predictive distributions of "sample path statistics" that might include, for example @@ -156,76 +159,141 @@ class AR1(NamedTuple): T1: int=100 -def AR1_simulate(ar1: AR1, T, subkey): +def AR1_simulate_past(ar1: AR1, key=key): """ - Simulate a realization of the AR(1) process for a given number of periods. + Simulate a realization of the AR(1) process for T0 periods. Parameters ---------- ar1 : AR1 - AR1 named tuple containing process parameters (rho, sigma, y0). - T : int - Number of time steps to simulate. - subkey : jax.random.PRNGKey + AR1 named tuple containing parameters (rho, sigma, y0, T0, T1). + key : jax.random.PRNGKey JAX random key for generating random noise. Returns ------- - y : jax.numpy.ndarray + initial_path : jax.numpy.ndarray Simulated path of the AR(1) process of length T. """ - # ...function code... - rho, sigma, y0 = ar.rho, ar1.sigma, ar1.y0 - # Allocate space and draw epsilons - y = jnp.empty(T) - eps = ar1.sigma * random.normal(subkey, (T,)) - - # Initial condition and step forward - y.at[0].set(y0) - for t in range(1, T): - y.at[t].set(ar1.rho * y[t-1] + eps[t]) - - return jnp.asarray(y) - - -def plot_path(ar1, initial_path, key, ax): - """Plot the initial path and the preceding predictive densities""" + rho, sigma, y0, T0 = ar1.rho, ar1.sigma, ar1.y0, ar1.T0 + # Draw epsilons + eps = sigma * random.normal(key, (T0 - 1,)) + # Set step function + def ar1_step(y_prev, t_rho_eps): + rho, eps_t = t_rho_eps + y_t = rho * y_prev + eps_t + return y_t, y_t + + # Scan over time steps + _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0-1, rho), eps)) + # Concatenate initial value + initial_path = jnp.concatenate([jnp.array([y0]), y_seq]) + + return initial_path + + +def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): + """ + Simulate a realization of the AR(1) process for T1 periods. + + Parameters + ---------- + ar1 : AR1 + AR1 named tuple containing parameters (rho, sigma, y0, T0, T1). + y_T0 : float + Value of the process at time T0. + N: int + Number of paths to simulate. + key : jax.random.PRNGKey + JAX random key for generating random noise. + + Returns + ------- + future_path : jax.numpy.ndarray + Simulated N paths of the AR(1) process of length T1. + """ + rho, sigma, T1 = ar1.rho, ar1.sigma, ar1.T1 + + def single_path_scan(y_T0, subkey): + eps = sigma * random.normal(subkey, (T1,)) + def ar1_step(y_prev, t_rho_eps): + rho, eps_t = t_rho_eps + y_t = rho * y_prev + eps_t + return y_t, y_t + _, y = lax.scan(ar1_step, y_T0, (jnp.full(T1, rho), eps)) + return y + + # Split key to generate different paths + subkeys = random.split(key, num=N) + # Simulate N paths + future_path = jax.vmap(single_path_scan, in_axes=(None, 0))(y_T0, subkeys) + + return future_path + + +def plot_path(ar1, initial_path, future_path, ax, key=key): + """ + Plot the initial observed AR(1) path and simulated future paths, + along with predictive confidence intervals. + + Parameters + ---------- + ar1 : AR1 + AR1 named tuple containing process parameters (rho, sigma, T0, T1). + initial_path : array-like + Simulated initial path of the AR(1) process of length T0. + future_path : array-like + Simulated future paths of the AR(1) process, shape (N, T1). + ax : matplotlib.axes.Axes + Matplotlib axes object to plot on. + key : jax.random.PRNGKey, optional + JAX random key for reproducible sampling. + + Plots + ----- + - Initial path (historical data) + - Multiple simulated future paths + - 90% and 95% predictive confidence intervals + - Expected future path + """ rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 # Compute moments and confidence intervals y_T0 = initial_path[-1] - center = jnp.array([rho**j * y_T0 for j in range(T1)]) - vars = jnp.array( - [sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)] - ) - + j = jnp.arange(T1) + center = rho**j * y_T0 + vars = sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) + # 95% CI y_upper_c95 = center + 1.96 * jnp.sqrt(vars) y_lower_c95 = center - 1.96 * jnp.sqrt(vars) - + # 90% CI y_upper_c90 = center + 1.65 * jnp.sqrt(vars) y_lower_c90 = center - 1.65 * jnp.sqrt(vars) # Plot - ax.plot(jnp.arange(-T0 + 1, 1), initial_path) - ax.set_xlim([-T0, T1]) + ax.plot(jnp.arange(-T0, 0), initial_path) ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1) + # Choose 10 future paths to plot + index = random.choice( + key, jnp.arange(future_path.shape[0]), (10,), replace=False + ) + for i in index: + ax.plot(jnp.arange(T1), future_path[i, :], color='grey', alpha=.5) - # Simulate 10 future paths - subkeys = random.split(key, num=10) - for i in range(10): - ar1_n = AR1(rho=rho, sigma=sigma, y0=y_T0) - y_future = AR1_simulate(ar1, T1, subkeys[i]) - ax.plot(jnp.arange(T1), y_future, color='grey', alpha=.5) - - # Plot 90% CI + # Plot 90% and 95% CI ax.fill_between( jnp.arange(T1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' ) ax.fill_between( jnp.arange(T1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' ) - ax.plot(jnp.arange(T1), center, color='red', alpha=.7, label='expectation') - ax.legend(fontsize=12) + ax.plot( + jnp.arange(T1), center, color='red', alpha=.7, label='expectation' + ) + ax.set_xlim([-T0, T1]) + ax.set_xlabel("time", fontsize=13) + ax.set_ylabel("y", fontsize=13) + ax.legend(fontsize=8) ``` ```{code-cell} ipython3 @@ -233,17 +301,17 @@ def plot_path(ar1, initial_path, key, ax): mystnb: figure: caption: | - Initial path and predictive densities + Initial and predictive future paths name: fig_path --- ar1 = AR1(rho=0.9, sigma=1, y0=10, T0=100, T1=100) # Simulate -initial_path = AR1_simulate(ar1, ar1.T0, key) - +initial_path = AR1_simulate_past(ar1) +future_path = AR1_simulate_future(ar1, initial_path[-1]) # Plot fig, ax = plt.subplots(1, 1) -plot_path(ar1, initial_path, key, ax) +plot_path(ar1, initial_path, future_path, ax) plt.show() ``` @@ -348,39 +416,45 @@ mystnb: Posterior distributions (rho, sigma) name: fig_post --- -def draw_from_posterior(sample, size=1e5, bins=20, dis_plot=1): +def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): """Draw a sample of size from the posterior distribution.""" - # Start with priors - rho = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable rho - sigma = numpyro.sample('sigma', dis.HalfNormal(sqrt(10))) - - # Define likelihood recursively - for t in range(1, len(sample)) - # Expectation of y_t - mu = rho * sample[t-1] - - # Likelihood of the actual realization. - numpyro.sample(f'y_{t}', dist.Normal(mu, sigma), obs=sample[t]) + def model(data): + # Start with priors + rho = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable rho + sigma = numpyro.sample('sigma', dist.HalfNormal(jnp.sqrt(10))) + # Define likelihood recursively + for t in range(1, len(data)): + # Expectation of y_t + mu = rho * data[t-1] + # Likelihood of the actual realization. + numpyro.sample(f'y_{t}', dist.Normal(mu, sigma), obs=data[t]) # Compute posterior distribution of parameters nuts_kernel = NUTS(model) + # Define MCMC class to compute the posteriors mcmc = MCMC( - nuts_kernek, - num_warmup=5000, - num_sample=size, + nuts_kernel, + num_warmup=5000, + num_samples=size, progress_bar=False) - mcmc.run(random.PRNGKey(0), sample=sample) - + # Run MCMC + mcmc.run(key, data=data) + # Get posterior samples post_sample = { 'rho': mcmc.get_samples()['rho'], 'sigma': mcmc.get_samples()['sigma'] } - + # Plot posterior distributions if dis_plot==1: - sns.displot( - post_sample, kde=True, stat="density", bins=bins, height=5, aspect=1.5 - ) - # TODO:Check plot and return + fig, ax = plt.subplots(1, 2, figsize=(12, 5)) + sns.histplot( + post_sample['rho'], kde=True, stat="density", bins=bins, ax=ax[0] + ) + ax[0].set_xlabel(r"$\rho$") + sns.histplot( + post_sample['sigma'], kde=True, stat="density", bins=bins, ax=ax[1] + ) + ax[1].set_xlabel(r"$\sigma$") return post_sample @@ -398,61 +472,75 @@ These statistics were originally defined as random variables with respect to $\o Also, these two kinds of definitions are equivalent because $\omega$ determins path statistics only through $\{Y_t\}$ ```{code-cell} ipython3 -# Define statistics as functions of y, the realized path values. -def next_recession(y): - n = y.shape[0] - 3 - z = jnp.zeros(n, dtype=int) +def compute_path_statistics(initial_path, future_path): + """Compute path statistics for the AR(1) process.""" + # Concatenate the last two elements of initial path to identify recession + y = jnp.concatenate([initial_path[-2:], future_path]) + n = y.shape[0] + def step(carry, i): + # identify recession + rec_cond = (y[i] < y[i-1]) & (y[i-1] < y[i-2]) + # identify severe recession + sev_cond = (y[i] - y[i-1] < -0.02) & (y[i-1] - y[i-2] < -0.02) + # identify positive turning point + up_cond = ( + ((y[i-2] > y[i-1]) & (y[i-1] > y[i]) & (y[i] < y[i+1]) & (y[i+1] < y[i+2])) + | ((y[i-1] > y[i]) & (y[i] > y[i+1]) & (y[i+1] < y[i+2]) & (y[i+2] < y[i+3])) + ) + # identify negative turning point + down_cond = ( + (y[i-2] < y[i-1]) & (y[i-1] < y[i]) & (y[i] > y[i+1]) & (y[i+1] > y[i+2]) + | ((y[i-1] < y[i]) & (y[i] < y[i+1]) & (y[i+1] > y[i+2]) & (y[i+2] > y[i+3])) + ) + # Convert to int + rec = jnp.where(rec_cond, 1, 0) + sev = jnp.where(sev_cond, 1, 0) + up = jnp.where(up_cond, 1, 0) + down = jnp.where(down_cond, 1, 0) + return carry, (rec, sev, up, down) - for i in range(n): - z.at[i].set(int(y[i] <= y[i+1] and y[i+1] > y[i+2] and y[i+2] > y[i+3])) - - if jnp.any(z) == False: - return 500 - else: - return jnp.where(z == 1)[0][0] + 1 - - -def next_severe_recession(y): - z = jnp.diff(y) - n = z.shape[0] + _, (rec_seq, sev_seq, up_seq, down_seq) = lax.scan(step, None, jnp.arange(2, n-3)) - sr = (z < -.02).astype(int) - indices = jnp.where(sr == 1)[0] + # Get the index of the first recession + next_recession = jnp.where( + jnp.any(rec_seq == 1), jnp.argmax(rec_seq == 1) + 1, len(y) + ) + next_severe_recession = jnp.where( + jnp.any(sev_seq == 1), jnp.argmax(sev_seq == 1) + 1, len(y) + ) + min_val_8q = jnp.min(future_path[:8]) # Minimum value in the next 8 periods + next_up_turn = jnp.where( + jnp.any(up_seq == 1), jnp.argmax(up_seq == 1) + 1, len(y) + ) + next_down_turn = jnp.where( + jnp.any(down_seq == 1), jnp.argmax(down_seq == 1) + 1, len(y) + ) - if len(indices) == 0: - return T1 - else: - return indices[0] + 1 + path_stats = (next_recession, next_severe_recession, min_val_8q, + next_up_turn, next_down_turn) + return path_stats -def minimum_value_8q(y): - return jnp.min(y[:8]) +def plot_path_stats(next_recession, next_severe_recession, min_val_8q, + next_up_turn, next_down_turn, ax): + """Plot the path statistics in subplots(3,2)""" + # ax[0, 0] is for paths of y + sns.histplot(next_recession, kde=True, stat='density', ax=ax[0, 1], alpha=.8) + ax[0, 1].set_xlabel("time until the next recession", fontsize=13) + sns.histplot( + next_severe_recession, kde=True, stat='density', ax=ax[1, 0], alpha=.8 + ) + ax[1, 0].set_xlabel("time until the next severe recession",fontsize=13) -def next_turning_point(y): - """ - Suppose that y is of length 6 - - y_{t-2}, y_{t-1}, y_{t}, y_{t+1}, y_{t+2}, y_{t+3} - - that is sufficient for determining the value of P/N - """ - - n = jnp.asarray(y).shape[0] - 4 - T = jnp.zeros(n, dtype=int) - - for i in range(n): - if ((y[i] > y[i+1]) and (y[i+1] > y[i+2]) and - (y[i+2] < y[i+3]) and (y[i+3] < y[i+4])): - T[i] = 1 - elif ((y[i] < y[i+1]) and (y[i+1] < y[i+2]) and - (y[i+2] > y[i+3]) and (y[i+3] > y[i+4])): - T[i] = -1 - - up_turn = jnp.where(T == 1)[0][0] + 1 if (1 in T) == True else T1 - down_turn = jnp.where(T == -1)[0][0] + 1 if (-1 in T) == True else T1 + sns.histplot(min_val_8q, kde=True, stat='density', ax=ax[1, 1], alpha=.8) + ax[1, 1].set_xlabel("minimum value in next 8 periods", fontsize=13) + + sns.histplot(next_up_turn, kde=True, stat='density', ax=ax[2, 0], alpha=.8) + ax[2, 0].set_xlabel("time until the next positive turn", fontsize=13) - return up_turn, down_turn + sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.8) + ax[2, 1].set_xlabel("time until the next negative turn", fontsize=13) ``` ## Original Wecker Method @@ -467,7 +555,7 @@ mystnb: Distributions of statistics by Wecker's method name: fig_wecker --- -def plot_Wecker(ar1: AR1, initial_path, N, ax): +def plot_Wecker(ar1: AR1, initial_path, ax, N=1000): """ Plot the predictive distributions from "pure" Wecker's method. @@ -480,98 +568,29 @@ def plot_Wecker(ar1: AR1, initial_path, N, ax): N : int Number of future sample paths to simulate for predictive distributions. """ - rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 - # Store outcomes - next_reces = jnp.zeros(N) - severe_rec = jnp.zeros(N) - min_vals = jnp.zeros(N) - next_up_turn, next_down_turn = jnp.zeros(N), jnp.zeros(N) - - # Compute confidence intervals + # Plot simulated initial and future paths y_T0 = initial_path[-1] - center = jnp.array([rho**j * yT0 for j in range(T1)]) - vars = jnp.array( - [sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)] - ) - y_upper_c95 = center + 1.96 * jnp.sqrt(vars) - y_lower_c95 = center - 1.96 * jnp.sqrt(vars) - - y_upper_c90 = center + 1.65 * jnp.sqrt(vars) - y_lower_c90 = center - 1.65 * jnp.sqrt(vars) - - # Plot - ax[0, 0].set_title("Initial path and predictive densities", fontsize=15) - ax[0, 0].plot(jnp.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 - ax[0, 0].fill_between(jnp.arange(T1), y_upper_c95, y_lower_c95, alpha=.3) - ax[0, 0].fill_between(jnp.arange(T1), y_upper_c90, y_lower_c90, alpha=.35) - ax[0, 0].plot(jnp.arange(T1), center, color='red', alpha=.7) - - # Simulate future paths - subkeys = random.split(key, num=N) - for n in range(N): - ar1_n = AR1(rho=rho, sigma=sigma, y0=y_T0) - sim_path = AR1_simulate(ar1_n, T1, subkeys[n]) - next_reces[n] = next_recession(jnp.hstack([initial_path[-3:-1], sim_path])) - severe_rec[n] = next_severe_recession(sim_path) - min_vals[n] = minimum_value_8q(sim_path) - next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) - - if n%(N/10) == 0: - ax[0, 0].plot(jnp.arange(T1), sim_path, color='gray', alpha=.3, lw=1) - - # 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 (time until the next recession)", - fontsize=13 - ) - - sns.histplot( - severe_rec, kde=False, stat='density', - ax=ax[1, 0], binwidth=0.9, alpha=.7, label='True parameters' - ) - ax[1, 0].set_title( - r"Predictive distribution (stopping time of growth$<-2\%$)", - fontsize=13 - ) - - sns.histplot( - min_vals, kde=True, stat='density', - ax=ax[1, 1], alpha=.8, label='True parameters' - ) - ax[1, 1].set_title( - "Predictive distribution (minimum value in next 8 periods)", - fontsize=13 - ) - - sns.histplot( - next_up_turn, kde=True, stat='density', - ax=ax[2, 0], alpha=.8, label='True parameters' - ) - ax[2, 0].set_title( - "Predictive distribution (time until the next positive turn)", - fontsize=13 - ) - - sns.histplot( - next_down_turn, kde=True, stat='density', - ax=ax[2, 1], alpha=.8, label='True parameters' - ) - ax[2, 1].set_title( - "Predictive distribution (time until the next negative turn)", - fontsize=13 - ) + future_path = AR1_simulate_future(ar1, y_T0, N=N) + plot_path(ar1, initial_path, future_path, ax[0, 0]) + + # Simulate future paths and compute statistics + def step(carry, n): + future_temp = future_path[n, :] + (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + ) = compute_path_statistics(initial_path, future_temp) + + return carry, (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn) + + _, (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + ) = lax.scan(step, None, jnp.arange(N)) + + # Plot path statistics + plot_path_stats(next_reces, severe_rec, min_vals, + next_up_turn, next_down_turn, ax) -fig, ax = plt.subplots(3, 2) -plot_Wecker(ar1, initial_path, 1000) +fig, ax = plt.subplots(3, 2, figsize=(15, 12)) +plot_Wecker(ar1, initial_path, ax) plt.show() ``` @@ -590,12 +609,14 @@ mystnb: Distributions of statistics by extended Wecker's method name: fig_extend_wecker --- -def plot_extended_Wecker(ar1, post_samples, initial_path, N, ax): - """ - Plot the extended Wecker's predictive distribution - """ - rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 - # Select a sample +def plot_extended_Wecker( + ar1: AR1, post_samples, initial_path, ax, N=1000 + ): + """Plot the extended Wecker's predictive distribution""" + y0, T1 = ar1.y0, ar1.T1 + y_T0 = initial_path[-1] + + # Select a parameter sample index = random.choice( key, jnp.arange(len(post_samples['rho'])), (N + 1,), replace=False ) @@ -603,87 +624,35 @@ def plot_extended_Wecker(ar1, post_samples, initial_path, N, ax): sigma_sample = post_samples['sigma'][index] # Store outcomes + future_path = jnp.zeros((N, T1)) next_reces = jnp.zeros(N) severe_rec = jnp.zeros(N) min_vals = jnp.zeros(N) next_up_turn, next_down_turn = jnp.zeros(N), jnp.zeros(N) - # Plot - ax[0, 0].set_title( - "Initial path and future paths simulated from posterior draws", - fontsize=15 - ) - ax[0, 0].plot(jnp.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 + # Compute path statistics subkeys = random.split(key, num=N) - for n in range(N): - ar1_n = AR1( - rho=rho_sample[n], sigma=sigma_sample[n], y0= initial_path[-1] - ) - sim_path = AR1_simulate(ar1_n, T1, subkeys[n]) - next_reces[n] = next_recession( - jnp.hstack([initial_path[-3:-1], sim_path]) - ) - severe_rec[n] = next_severe_recession(sim_path) - min_vals[n] = minimum_value_8q(sim_path) - next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path) - - if n % (N / 10) == 0: - ax[0, 0].plot(jnp.arange(T1), sim_path, color='gray', alpha=.3, lw=1) - - # 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 (time until the next recession)", - fontsize=13) - - sns.histplot( - severe_rec, kde=False, stat='density', - ax=ax[1, 0], binwidth=.9, alpha=.6, color=colors[1], - label='Sampling from posterior' - ) - ax[1, 0].set_title( - r"Predictive distribution (stopping time of growth$<-2\%$)", - fontsize=13 - ) - - sns.histplot( - min_vals, kde=True, stat='density', - ax=ax[1, 1], alpha=.6, color=colors[1], - label='Sampling from posterior' - ) - ax[1, 1].set_title( - "Predictive distribution (minimum value in the next 8 periods)", - fontsize=13) - - sns.histplot( - next_up_turn, kde=True, stat='density', - ax=ax[2, 0], alpha=.6, color=colors[1], - label='Sampling from posterior') - ax[2, 0].set_title( - "Predictive distribution of time until the next positive turn", - fontsize=13 - ) - - 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) -plot_extended_Wecker(ar1, post_samples, initial_path, 1000, ax) + def step(carry, n): + ar1_n = AR1(rho=rho_sample[n], sigma=sigma_sample[n], y0=y0) + future_temp = AR1_simulate_future( + ar1_n, y_T0, N=1, key=subkeys[n] + ).reshape(-1) + (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + ) = compute_path_statistics(initial_path, future_temp) + return carry, (future_temp, next_reces, severe_rec, + min_vals, next_up_turn, next_down_turn) + + _, (future_path, next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + ) = jax.lax.scan(step, None, jnp.arange(N)) + + # Plot simulated initial and future paths + plot_path(ar1, initial_path, future_path, ax[0, 0]) + # Plot path statistics + plot_path_stats(next_reces, severe_rec, min_vals, + next_up_turn, next_down_turn, ax) + +fig, ax = plt.subplots(3, 2, figsize=(12, 15)) +plot_extended_Wecker(ar1, post_samples, initial_path, ax) plt.show() ``` @@ -699,10 +668,10 @@ mystnb: Comparison between two methods name: fig_wecker --- -fig, ax = plt.subplots(3, 2) -plot_Wecker(ar1, initial_path, 1000, ax) +fig, ax = plt.subplots(3, 2, figsize=(12, 15)) +plot_Wecker(ar1, initial_path, ax) ax[0, 0].clear() -plot_extended_Wecker(ar1, post_samples, initial_path, 1000, ax) +plot_extended_Wecker(ar1, post_samples, initial_path, ax) plt.legend() plt.show() ``` From cfbc63209610341d3a22861166bf2188f72b16ee Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Mon, 1 Sep 2025 10:50:44 +0800 Subject: [PATCH 05/16] fix typos and add explanations --- lectures/ar1_turningpts.md | 104 ++++++++++++++++++++++--------------- 1 file changed, 62 insertions(+), 42 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 2dac59d9a..ba2f6443e 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -40,7 +40,7 @@ We consider two sorts of statistics: To investigate sample path properties we'll use a simulation procedure recommended by Wecker {cite}`wecker1979predicting`. -To acknowledge uncertainty about parameters, we'll deploy `pymc` to construct a Bayesian joint posterior distribution for unknown parameters. +To acknowledge uncertainty about parameters, we'll deploy `numpyro` to construct a Bayesian joint posterior distribution for unknown parameters. Let's start with some imports. @@ -110,7 +110,7 @@ Predictive distribution {eq}`ar1-tp-eq3` assumes that parameters $(\rho,\sigma)$ Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows that $\{y_t\}$ is a AR(1) process when $(\rho, \sigma)$ are given. -We also want to compute some predictive distributions of "sample path statistics" that might include, for example +We also want to compute some predictive distributions of "sample path statistics" that might include, for example - the time until the next "recession", - the minimum value of $Y$ over the next 8 periods, @@ -152,8 +152,8 @@ class AR1(NamedTuple): T1 : int, optional Length of the future path to simulate (default is 100). """ - rho:float - sigma:float + rho: float + sigma: float y0: float T0: int=100 T1: int=100 @@ -173,11 +173,11 @@ def AR1_simulate_past(ar1: AR1, key=key): Returns ------- initial_path : jax.numpy.ndarray - Simulated path of the AR(1) process of length T. + Simulated path of the AR(1) process and the initial y0. """ rho, sigma, y0, T0 = ar1.rho, ar1.sigma, ar1.y0, ar1.T0 # Draw epsilons - eps = sigma * random.normal(key, (T0 - 1,)) + eps = sigma * random.normal(key, (T0,)) # Set step function def ar1_step(y_prev, t_rho_eps): rho, eps_t = t_rho_eps @@ -185,7 +185,7 @@ def AR1_simulate_past(ar1: AR1, key=key): return y_t, y_t # Scan over time steps - _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0-1, rho), eps)) + _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0, rho), eps)) # Concatenate initial value initial_path = jnp.concatenate([jnp.array([y0]), y_seq]) @@ -241,7 +241,7 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): ar1 : AR1 AR1 named tuple containing process parameters (rho, sigma, T0, T1). initial_path : array-like - Simulated initial path of the AR(1) process of length T0. + Simulated initial path of the AR(1) process, shape(T0+1,). future_path : array-like Simulated future paths of the AR(1) process, shape (N, T1). ax : matplotlib.axes.Axes @@ -260,7 +260,7 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): # Compute moments and confidence intervals y_T0 = initial_path[-1] - j = jnp.arange(T1) + j = jnp.arange(1, T1+1) center = rho**j * y_T0 vars = sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) # 95% CI @@ -271,24 +271,24 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): y_lower_c90 = center - 1.65 * jnp.sqrt(vars) # Plot - ax.plot(jnp.arange(-T0, 0), initial_path) + ax.plot(jnp.arange(-T0, 1), initial_path) ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1) - # Choose 10 future paths to plot + # Chooise 10 future paths to plot index = random.choice( key, jnp.arange(future_path.shape[0]), (10,), replace=False ) for i in index: - ax.plot(jnp.arange(T1), future_path[i, :], color='grey', alpha=.5) + ax.plot(jnp.arange(1, T1+1), future_path[i, :], color='grey', alpha=.5) # Plot 90% and 95% CI ax.fill_between( - jnp.arange(T1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' + jnp.arange(1, T1+1), y_upper_c95, y_lower_c95, alpha=.3, label='95% CI' ) ax.fill_between( - jnp.arange(T1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' + jnp.arange(1, T1+1), y_upper_c90, y_lower_c90, alpha=.35, label='90% CI' ) ax.plot( - jnp.arange(T1), center, color='red', alpha=.7, label='expectation' + jnp.arange(1, T1+1), center, color='red', alpha=.7, label='expectation' ) ax.set_xlim([-T0, T1]) ax.set_xlabel("time", fontsize=13) @@ -325,13 +325,15 @@ Wecker {cite}`wecker1979predicting` proposed using simulation techniques to char He called these functions "path properties" to contrast them with properties of single data points. -He studied two special prospective path properties of a given series $\{y_t\}$. +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$. +* he defined a **"turning point"** to be the date of the second of two successive declines in $y$. -To examine this statistic, let $Z$ be an indicator process +For example, if $y_t(\omega)< y_{t-1}(\omega)< y_{t-2}(\omega)$, then $t$, then period $t$ is a turning point. + +To examine the *time until the next turning point*, let $Z$ be an indicator process $$ Z_t(\omega) := @@ -343,12 +345,18 @@ $$ Here $\omega \in Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow R$ gives $y_t$ according to $\omega$ and the AR(1) process. -Then the random variable **time until the next turning point** is defined as the following *stopping time* with respect to $Z$: +By Wecker's definition, period $t$ is a turning point, and $Y_{t-2}(\omega) \geq Y_{t-3}(\omega)$ excludes that period $t-1$ is a turning point. + +Then the random variable **time until the next turning point** is defined as the following *stopping time* with respect to $Z$: $$ W_t(\omega):= \inf \{ k\geq 1 \mid Z_{t+k}(\omega) = 1\} $$ +In the following code, we name this statistic *time until the next recession* to distinguish it from another concept of *turning point*. + +Moreover, the statistic *time until the next severe recession* is defined in a similar way, except the decline between periods is greater than $0.02$. + Wecker {cite}`wecker1979predicting` also studied **the minimum value of $Y$ over the next 8 quarters** which can be defined as the random variable. @@ -383,8 +391,12 @@ This is designed to express the event - "after one or two decrease(s), $Y$ will grow for two consecutive quarters" +The **negative turning point today or tomorrow** $N_t$ is defined in the same way. + Following {cite}`wecker1979predicting`, we can use simulations to calculate probabilities of $P_t$ and $N_t$ for each period $t$. +However, in the following code, we only use $T_{t+1}(\omega)=1$ to determine $P_t(\omega)$ and $N_t(\omega)$, because we only want to find out the first positive turning point. + ## A Wecker-Like Algorithm The procedure consists of the following steps: @@ -404,7 +416,7 @@ $$ ## Using Simulations to Approximate a Posterior Distribution -The next code cells use `pymc` to compute the time $t$ posterior distribution of $\rho, \sigma$. +The next code cells use `numpyro` to compute the time $t$ posterior distribution of $\rho, \sigma$. Note that in defining the likelihood function, we choose to condition on the initial value $y_0$. @@ -445,7 +457,7 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): 'sigma': mcmc.get_samples()['sigma'] } # Plot posterior distributions - if dis_plot==1: + if dis_plot == 1: fig, ax = plt.subplots(1, 2, figsize=(12, 5)) sns.histplot( post_sample['rho'], kde=True, stat="density", bins=bins, ax=ax[0] @@ -469,28 +481,30 @@ Our next step is to prepare Python code to compute our sample path statistics. These statistics were originally defined as random variables with respect to $\omega$, but here we use $\{Y_t\}$ as the argument because $\omega$ is implicit. -Also, these two kinds of definitions are equivalent because $\omega$ determins path statistics only through $\{Y_t\}$ +These two kinds of definitions are equivalent because $\omega$ determins path statistics only through $\{Y_t\}$. + +Moreover, we ignore all equality in the definitions, as equality occurs with zero probablity for countinuous random variables. ```{code-cell} ipython3 def compute_path_statistics(initial_path, future_path): """Compute path statistics for the AR(1) process.""" # Concatenate the last two elements of initial path to identify recession - y = jnp.concatenate([initial_path[-2:], future_path]) + y = jnp.concatenate([initial_path[-3:], future_path]) n = y.shape[0] def step(carry, i): # identify recession - rec_cond = (y[i] < y[i-1]) & (y[i-1] < y[i-2]) + rec_cond = (y[i] < y[i-1]) & (y[i-1] < y[i-2]) & (y[i-2] > y[i-3]) # identify severe recession - sev_cond = (y[i] - y[i-1] < -0.02) & (y[i-1] - y[i-2] < -0.02) - # identify positive turning point + sev_cond = ( + (y[i] - y[i-1] < -0.02) & (y[i-1] - y[i-2] < -0.02) & (y[i-2] > y[i-3]) + ) + # identify positive turning point. up_cond = ( - ((y[i-2] > y[i-1]) & (y[i-1] > y[i]) & (y[i] < y[i+1]) & (y[i+1] < y[i+2])) - | ((y[i-1] > y[i]) & (y[i] > y[i+1]) & (y[i+1] < y[i+2]) & (y[i+2] < y[i+3])) + (y[i-2] > y[i-1]) & (y[i-1] > y[i]) & (y[i] < y[i+1]) & (y[i+1] < y[i+2]) ) # identify negative turning point down_cond = ( (y[i-2] < y[i-1]) & (y[i-1] < y[i]) & (y[i] > y[i+1]) & (y[i+1] > y[i+2]) - | ((y[i-1] < y[i]) & (y[i] < y[i+1]) & (y[i+1] > y[i+2]) & (y[i+2] > y[i+3])) ) # Convert to int rec = jnp.where(rec_cond, 1, 0) @@ -499,30 +513,37 @@ def compute_path_statistics(initial_path, future_path): down = jnp.where(down_cond, 1, 0) return carry, (rec, sev, up, down) - _, (rec_seq, sev_seq, up_seq, down_seq) = lax.scan(step, None, jnp.arange(2, n-3)) + _, (rec_seq, sev_seq, up_seq, down_seq) = lax.scan(step, None, jnp.arange(3, n-3)) - # Get the index of the first recession + # Get the time until the first recession next_recession = jnp.where( jnp.any(rec_seq == 1), jnp.argmax(rec_seq == 1) + 1, len(y) ) next_severe_recession = jnp.where( jnp.any(sev_seq == 1), jnp.argmax(sev_seq == 1) + 1, len(y) ) - min_val_8q = jnp.min(future_path[:8]) # Minimum value in the next 8 periods + # Minimum value in the next 8 periods + min_val_8q = jnp.min(future_path[:8]) + # Get the time until the first turning point. next_up_turn = jnp.where( - jnp.any(up_seq == 1), jnp.argmax(up_seq == 1) + 1, len(y) + jnp.any(up_seq == 1), + jnp.maximum(jnp.argmax(up_seq == 1), 1), # Exclude 0 return + len(y) ) next_down_turn = jnp.where( - jnp.any(down_seq == 1), jnp.argmax(down_seq == 1) + 1, len(y) + jnp.any(down_seq == 1), + jnp.maximum(jnp.argmax(down_seq == 1), 1), + len(y) ) - - path_stats = (next_recession, next_severe_recession, min_val_8q, - next_up_turn, next_down_turn) + path_stats = ( + next_recession, next_severe_recession, min_val_8q, + next_up_turn, next_down_turn + ) return path_stats def plot_path_stats(next_recession, next_severe_recession, min_val_8q, - next_up_turn, next_down_turn, ax): + next_up_turn, next_down_turn, ax): """Plot the path statistics in subplots(3,2)""" # ax[0, 0] is for paths of y sns.histplot(next_recession, kde=True, stat='density', ax=ax[0, 1], alpha=.8) @@ -531,7 +552,7 @@ def plot_path_stats(next_recession, next_severe_recession, min_val_8q, sns.histplot( next_severe_recession, kde=True, stat='density', ax=ax[1, 0], alpha=.8 ) - ax[1, 0].set_xlabel("time until the next severe recession",fontsize=13) + ax[1, 0].set_xlabel("time until the next severe recession", fontsize=13) sns.histplot(min_val_8q, kde=True, stat='density', ax=ax[1, 1], alpha=.8) ax[1, 1].set_xlabel("minimum value in next 8 periods", fontsize=13) @@ -666,12 +687,11 @@ mystnb: figure: caption: | Comparison between two methods - name: fig_wecker + name: fig_compare_wecker --- fig, ax = plt.subplots(3, 2, figsize=(12, 15)) plot_Wecker(ar1, initial_path, ax) ax[0, 0].clear() plot_extended_Wecker(ar1, post_samples, initial_path, ax) -plt.legend() plt.show() ``` From e09504a606d7e20a689a7303baf73e2bd50ea525 Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Mon, 1 Sep 2025 14:28:48 +0800 Subject: [PATCH 06/16] fix typo --- lectures/ar1_turningpts.md | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index ba2f6443e..f46435a2e 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -343,7 +343,7 @@ Z_t(\omega) := \end{cases} $$ -Here $\omega \in Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow R$ gives $y_t$ according to $\omega$ and the AR(1) process. +Here $\omega \in \Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow R$ gives $y_t$ according to $\omega$ and the AR(1) process. By Wecker's definition, period $t$ is a turning point, and $Y_{t-2}(\omega) \geq Y_{t-3}(\omega)$ excludes that period $t-1$ is a turning point. @@ -353,7 +353,7 @@ $$ W_t(\omega):= \inf \{ k\geq 1 \mid Z_{t+k}(\omega) = 1\} $$ -In the following code, we name this statistic *time until the next recession* to distinguish it from another concept of *turning point*. +In the following code, we name this statistic as *time until the next recession* to distinguish it from another concept of *turning point*. Moreover, the statistic *time until the next severe recession* is defined in a similar way, except the decline between periods is greater than $0.02$. @@ -644,13 +644,6 @@ def plot_extended_Wecker( rho_sample = post_samples['rho'][index] sigma_sample = post_samples['sigma'][index] - # Store outcomes - future_path = jnp.zeros((N, T1)) - next_reces = jnp.zeros(N) - severe_rec = jnp.zeros(N) - min_vals = jnp.zeros(N) - next_up_turn, next_down_turn = jnp.zeros(N), jnp.zeros(N) - # Compute path statistics subkeys = random.split(key, num=N) def step(carry, n): From d06a556e6148200acaeb40e3cfe0d06b526055e8 Mon Sep 17 00:00:00 2001 From: xuanguang-li Date: Tue, 2 Sep 2025 17:46:53 +0800 Subject: [PATCH 07/16] check grammar and style --- lectures/ar1_turningpts.md | 112 +++++++++++++++++++++++-------------- 1 file changed, 71 insertions(+), 41 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index f46435a2e..76036b8ae 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -22,11 +22,11 @@ kernelspec: !pip install numpyro jax ``` -This lecture describes methods for forecasting statistics that are functions of future values of a univariate autogressive process. +This lecture describes methods for forecasting statistics that are functions of future values of a univariate autoregressive process. The methods are designed to take into account two possible sources of uncertainty about these statistics: -- random shocks that impinge of the transition law +- random shocks that impinge on the transition law - uncertainty about the parameter values of the AR(1) process @@ -90,7 +90,7 @@ $$ f(y_{t+j} | y_{t}; \rho, \sigma) \sim {\mathcal N}\left(\rho^j y_{t}, \sigma^2 \frac{1 - \rho^{2j}}{1 - \rho^2} \right) $$ (ar1-tp-eq3) -The predictive distribution {eq}`ar1-tp-eq3` that assumes that the parameters $\rho, \sigma$ are known, which we express +The predictive distribution {eq}`ar1-tp-eq3` assumes that the parameters $\rho, \sigma$ are known, which we express by conditioning on them. We also want to compute a predictive distribution that does not condition on $\rho,\sigma$ but instead takes account of our uncertainty about them. @@ -108,7 +108,7 @@ $$ (ar1-tp-eq4) Predictive distribution {eq}`ar1-tp-eq3` assumes that parameters $(\rho,\sigma)$ are known. -Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows that $\{y_t\}$ is a AR(1) process when $(\rho, \sigma)$ are given. +Predictive distribution {eq}`ar1-tp-eq4` assumes that parameters $(\rho,\sigma)$ are uncertain, but have a known probability distribution $\pi_t(\rho,\sigma | y^t )$. Notice the second equality follows because $\{y_t\}$ is an AR(1) process when $(\rho, \sigma)$ are given. We also want to compute some predictive distributions of "sample path statistics" that might include, for example @@ -121,18 +121,18 @@ To accomplish that for situations in which we are uncertain about parameter valu - 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; -- for each draw $n=0,1,...,N$, simulate a "future path" of length $T_1$ with parameters $\left(\rho_n,\sigma_n\right)$ and compute our three "sample path statistics"; +- for each draw $n=0,1,...,N$, simulate a "future path" of length $T_1$ with parameters $\left(\rho_n,\sigma_n\right)$ and compute our "sample path statistics"; - finally, plot the desired statistics from the $N$ samples as an empirical distribution. ## Implementation -First, we'll simulate a sample path from which to launch our forecasts. +First, we'll simulate a sample path from which to launch our forecasts. In addition to plotting the sample path, under the assumption that the true parameter values are known, -we'll plot $.9$ and $.95$ coverage intervals using conditional distribution +we'll plot $0.9$ and $0.95$ coverage intervals using conditional distribution {eq}`ar1-tp-eq3` described above. -We'll also plot a bunch of samples of sequences of future values and watch where they fall relative to the coverage interval. +We'll also plot a bunch of samples of sequences of future values and watch where they fall relative to the coverage interval. ```{code-cell} ipython3 class AR1(NamedTuple): @@ -178,6 +178,7 @@ def AR1_simulate_past(ar1: AR1, key=key): rho, sigma, y0, T0 = ar1.rho, ar1.sigma, ar1.y0, ar1.T0 # Draw epsilons eps = sigma * random.normal(key, (T0,)) + # Set step function def ar1_step(y_prev, t_rho_eps): rho, eps_t = t_rho_eps @@ -186,12 +187,16 @@ def AR1_simulate_past(ar1: AR1, key=key): # Scan over time steps _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0, rho), eps)) + # Concatenate initial value initial_path = jnp.concatenate([jnp.array([y0]), y_seq]) return initial_path +``` +Now we define the simulation function that generates a realization of the AR(1) process for future $T1$ periods. +```{code-cell} ipython3 def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): """ Simulate a realization of the AR(1) process for T1 periods. @@ -216,6 +221,7 @@ def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): def single_path_scan(y_T0, subkey): eps = sigma * random.normal(subkey, (T1,)) + def ar1_step(y_prev, t_rho_eps): rho, eps_t = t_rho_eps y_t = rho * y_prev + eps_t @@ -225,12 +231,16 @@ def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): # Split key to generate different paths subkeys = random.split(key, num=N) + # Simulate N paths future_path = jax.vmap(single_path_scan, in_axes=(None, 0))(y_T0, subkeys) return future_path +``` +The following function plots the initial observed AR(1) path and simulated future paths along with predictive confidence intervals. +```{code-cell} ipython3 def plot_path(ar1, initial_path, future_path, ax, key=key): """ Plot the initial observed AR(1) path and simulated future paths, @@ -263,9 +273,11 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): j = jnp.arange(1, T1+1) center = rho**j * y_T0 vars = sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) + # 95% CI y_upper_c95 = center + 1.96 * jnp.sqrt(vars) y_lower_c95 = center - 1.96 * jnp.sqrt(vars) + # 90% CI y_upper_c90 = center + 1.65 * jnp.sqrt(vars) y_lower_c90 = center - 1.65 * jnp.sqrt(vars) @@ -273,7 +285,8 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): # Plot ax.plot(jnp.arange(-T0, 1), initial_path) ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1) - # Chooise 10 future paths to plot + + # Choose 10 future paths to plot index = random.choice( key, jnp.arange(future_path.shape[0]), (10,), replace=False ) @@ -309,6 +322,7 @@ ar1 = AR1(rho=0.9, sigma=1, y0=10, T0=100, T1=100) # Simulate initial_path = AR1_simulate_past(ar1) future_path = AR1_simulate_future(ar1, initial_path[-1]) + # Plot fig, ax = plt.subplots(1, 1) plot_path(ar1, initial_path, future_path, ax) @@ -321,17 +335,17 @@ https://python.quantecon.org/perm_income_cons.html ## Predictive Distributions of Path Properties -Wecker {cite}`wecker1979predicting` proposed using simulation techniques to characterize predictive distribution of some statistics that are non-linear functions of $y$. +Wecker {cite}`wecker1979predicting` proposed using simulation techniques to characterize the predictive distribution of some statistics that are non-linear functions of $y$. He called these functions "path properties" to contrast them with properties of single data points. -He studied two special prospective path properties of a given series $\{y_t\}$. +He studied two special prospective path properties of a given series $\{y_t\}$. 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$. -For example, if $y_t(\omega)< y_{t-1}(\omega)< y_{t-2}(\omega)$, then $t$, then period $t$ is a turning point. +For example, if $y_t(\omega)< y_{t-1}(\omega)< y_{t-2}(\omega)$, then period $t$ is a turning point. To examine the *time until the next turning point*, let $Z$ be an indicator process @@ -343,9 +357,9 @@ Z_t(\omega) := \end{cases} $$ -Here $\omega \in \Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow R$ gives $y_t$ according to $\omega$ and the AR(1) process. +Here $\omega \in \Omega$ is a sequence of events, and $Y_t: \Omega \rightarrow \mathbb{R}$ gives $y_t$ according to $\omega$ and the AR(1) process. -By Wecker's definition, period $t$ is a turning point, and $Y_{t-2}(\omega) \geq Y_{t-3}(\omega)$ excludes that period $t-1$ is a turning point. +By Wecker's definition, period $t$ is a turning point, and $Y_{t-2}(\omega) \geq Y_{t-3}(\omega)$ excludes the possibility that period $t-1$ is a turning point. Then the random variable **time until the next turning point** is defined as the following *stopping time* with respect to $Z$: @@ -353,12 +367,12 @@ $$ W_t(\omega):= \inf \{ k\geq 1 \mid Z_{t+k}(\omega) = 1\} $$ -In the following code, we name this statistic as *time until the next recession* to distinguish it from another concept of *turning point*. +In the following code, we name this statistic *time until the next recession* to distinguish it from another concept of *turning point*. Moreover, the statistic *time until the next severe recession* is defined in a similar way, except the decline between periods is greater than $0.02$. -Wecker {cite}`wecker1979predicting` also studied **the minimum value of $Y$ over the next 8 quarters** -which can be defined as the random variable. +Wecker {cite}`wecker1979predicting` also studied **the minimum value of $Y$ over the next 8 quarters**, +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)\} @@ -389,29 +403,29 @@ $$ This is designed to express the event -- "after one or two decrease(s), $Y$ will grow for two consecutive quarters" +- "after one or two decreases, $Y$ will grow for two consecutive quarters" The **negative turning point today or tomorrow** $N_t$ is defined in the same way. -Following {cite}`wecker1979predicting`, we can use simulations to calculate probabilities of $P_t$ and $N_t$ for each period $t$. +Following {cite}`wecker1979predicting`, we can use simulations to calculate probabilities of $P_t$ and $N_t$ for each period $t$. -However, in the following code, we only use $T_{t+1}(\omega)=1$ to determine $P_t(\omega)$ and $N_t(\omega)$, because we only want to find out the first positive turning point. +However, in the following code, we only use $T_{t+1}(\omega)=1$ to determine $P_t(\omega)$ and $N_t(\omega)$, because we only want to find the first positive turning point. ## A Wecker-Like Algorithm The procedure consists of the following steps: -* index a sample path by $\omega_i$ +- index a sample path by $\omega_i$ -* from a given date $t$, simulate $I$ sample paths of length $N$ +- from a given date $t$, simulate $I$ sample paths of length $N$ $$ Y(\omega_i) = \left\{ Y_{t+1}(\omega_i), Y_{t+2}(\omega_i), \dots, Y_{t+N}(\omega_i)\right\}_{i=1}^I $$ -* for each path $\omega_i$, compute the associated value of $W_t(\omega_i), W_{t+1}(\omega_i), \dots , W_{t+N}$ +- for each path $\omega_i$, compute the associated value of $W_t(\omega_i), W_{t+1}(\omega_i), \dots , W_{t+N}$ -* consider the sets $\{W_t(\omega_i)\}^{I}_{i=1}, \ \{W_{t+1}(\omega_i)\}^{I}_{i=1}, \ \dots, \ \{W_{t+N}(\omega_i)\}^{I}_{i=1}$ as samples from the predictive distributions $f(W_{t+1} \mid y_t, y_{t-1}, \dots , y_0)$, $f(W_{t+2} \mid y_t, y_{t-1}, \dots , y_0)$, $\dots$, $f(W_{t+N} \mid y_t, y_{t-1}, \dots , y_0)$. +- consider the sets $\{W_t(\omega_i)\}^{I}_{i=1}, \ \{W_{t+1}(\omega_i)\}^{I}_{i=1}, \ \dots, \ \{W_{t+N}(\omega_i)\}^{I}_{i=1}$ as samples from the predictive distributions $f(W_{t+1} \mid y_t, y_{t-1}, \dots , y_0)$, $f(W_{t+2} \mid y_t, y_{t-1}, \dots , y_0)$, $\dots$, $f(W_{t+N} \mid y_t, y_{t-1}, \dots , y_0)$. ## Using Simulations to Approximate a Posterior Distribution @@ -434,28 +448,34 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): # Start with priors rho = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable rho sigma = numpyro.sample('sigma', dist.HalfNormal(jnp.sqrt(10))) + # Define likelihood recursively for t in range(1, len(data)): # Expectation of y_t mu = rho * data[t-1] + # Likelihood of the actual realization. numpyro.sample(f'y_{t}', dist.Normal(mu, sigma), obs=data[t]) # Compute posterior distribution of parameters nuts_kernel = NUTS(model) + # Define MCMC class to compute the posteriors mcmc = MCMC( nuts_kernel, num_warmup=5000, num_samples=size, progress_bar=False) + # Run MCMC mcmc.run(key, data=data) + # Get posterior samples post_sample = { 'rho': mcmc.get_samples()['rho'], 'sigma': mcmc.get_samples()['sigma'] } + # Plot posterior distributions if dis_plot == 1: fig, ax = plt.subplots(1, 2, figsize=(12, 5)) @@ -481,9 +501,9 @@ Our next step is to prepare Python code to compute our sample path statistics. These statistics were originally defined as random variables with respect to $\omega$, but here we use $\{Y_t\}$ as the argument because $\omega$ is implicit. -These two kinds of definitions are equivalent because $\omega$ determins path statistics only through $\{Y_t\}$. +These two kinds of definitions are equivalent because $\omega$ determines path statistics only through $\{Y_t\}$. -Moreover, we ignore all equality in the definitions, as equality occurs with zero probablity for countinuous random variables. +Moreover, we ignore all equality in the definitions, as equality occurs with zero probability for continuous random variables. ```{code-cell} ipython3 def compute_path_statistics(initial_path, future_path): @@ -491,6 +511,7 @@ def compute_path_statistics(initial_path, future_path): # Concatenate the last two elements of initial path to identify recession y = jnp.concatenate([initial_path[-3:], future_path]) n = y.shape[0] + def step(carry, i): # identify recession rec_cond = (y[i] < y[i-1]) & (y[i-1] < y[i-2]) & (y[i-2] > y[i-3]) @@ -498,14 +519,17 @@ def compute_path_statistics(initial_path, future_path): sev_cond = ( (y[i] - y[i-1] < -0.02) & (y[i-1] - y[i-2] < -0.02) & (y[i-2] > y[i-3]) ) - # identify positive turning point. + + # identify positive turning point up_cond = ( (y[i-2] > y[i-1]) & (y[i-1] > y[i]) & (y[i] < y[i+1]) & (y[i+1] < y[i+2]) ) + # identify negative turning point down_cond = ( (y[i-2] < y[i-1]) & (y[i-1] < y[i]) & (y[i] > y[i+1]) & (y[i+1] > y[i+2]) ) + # Convert to int rec = jnp.where(rec_cond, 1, 0) sev = jnp.where(sev_cond, 1, 0) @@ -513,7 +537,7 @@ def compute_path_statistics(initial_path, future_path): down = jnp.where(down_cond, 1, 0) return carry, (rec, sev, up, down) - _, (rec_seq, sev_seq, up_seq, down_seq) = lax.scan(step, None, jnp.arange(3, n-3)) + _, (rec_seq, sev_seq, up_seq, down_seq) = lax.scan(step, None, jnp.arange(3, n-2)) # Get the time until the first recession next_recession = jnp.where( @@ -522,9 +546,11 @@ def compute_path_statistics(initial_path, future_path): next_severe_recession = jnp.where( jnp.any(sev_seq == 1), jnp.argmax(sev_seq == 1) + 1, len(y) ) + # Minimum value in the next 8 periods min_val_8q = jnp.min(future_path[:8]) - # Get the time until the first turning point. + + # Get the time until the first turning point next_up_turn = jnp.where( jnp.any(up_seq == 1), jnp.maximum(jnp.argmax(up_seq == 1), 1), # Exclude 0 return @@ -540,8 +566,10 @@ def compute_path_statistics(initial_path, future_path): next_up_turn, next_down_turn ) return path_stats +``` +The following function creates visualizations of the path statistics in a subplot grid. - +```{code-cell} ipython3 def plot_path_stats(next_recession, next_severe_recession, min_val_8q, next_up_turn, next_down_turn, ax): """Plot the path statistics in subplots(3,2)""" @@ -597,16 +625,16 @@ def plot_Wecker(ar1: AR1, initial_path, ax, N=1000): # Simulate future paths and compute statistics def step(carry, n): future_temp = future_path[n, :] - (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn ) = compute_path_statistics(initial_path, future_temp) - return carry, (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn) + return carry, (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn) - _, (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + _, (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn ) = lax.scan(step, None, jnp.arange(N)) # Plot path statistics - plot_path_stats(next_reces, severe_rec, min_vals, + plot_path_stats(next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn, ax) @@ -617,10 +645,10 @@ plt.show() ## Extended Wecker Method -Now we apply we apply our "extended" Wecker method based on predictive densities of $y$ defined by +Now 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 repeatedly 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 integration 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 --- @@ -646,23 +674,25 @@ def plot_extended_Wecker( # Compute path statistics subkeys = random.split(key, num=N) + def step(carry, n): ar1_n = AR1(rho=rho_sample[n], sigma=sigma_sample[n], y0=y0) future_temp = AR1_simulate_future( ar1_n, y_T0, N=1, key=subkeys[n] ).reshape(-1) - (next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn ) = compute_path_statistics(initial_path, future_temp) return carry, (future_temp, next_reces, severe_rec, - min_vals, next_up_turn, next_down_turn) + min_val_8q, next_up_turn, next_down_turn) - _, (future_path, next_reces, severe_rec, min_vals, next_up_turn, next_down_turn + _, (future_path, next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn ) = jax.lax.scan(step, None, jnp.arange(N)) # Plot simulated initial and future paths plot_path(ar1, initial_path, future_path, ax[0, 0]) + # Plot path statistics - plot_path_stats(next_reces, severe_rec, min_vals, + plot_path_stats(next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn, ax) fig, ax = plt.subplots(3, 2, figsize=(12, 15)) @@ -672,7 +702,7 @@ 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 differences 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 --- From dbd3605a2993ca4388198e52e8fc086eba16370f Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Thu, 4 Sep 2025 18:25:23 +0800 Subject: [PATCH 08/16] use Greek letters and arviz plot --- lectures/ar1_turningpts.md | 121 +++++++++++++++++++++---------------- 1 file changed, 69 insertions(+), 52 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 76036b8ae..af65a0862 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.16.4 + jupytext_version: 1.17.3 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -19,7 +19,7 @@ kernelspec: ```{code-cell} ipython3 :tags: [hide-output] -!pip install numpyro jax +!pip install numpyro jax arviz ``` This lecture describes methods for forecasting statistics that are functions of future values of a univariate autoregressive process. @@ -60,6 +60,9 @@ import numpyro import numpyro.distributions as dist from numpyro.infer import MCMC, NUTS +# arviz +import arviz as az + sns.set_style('white') colors = sns.color_palette() key = random.PRNGKey(0) @@ -141,9 +144,9 @@ class AR1(NamedTuple): Parameters ---------- - rho : float - Autoregressive coefficient, must satisfy |rho| < 1 for stationarity. - sigma : float + ρ : float + Autoregressive coefficient, must satisfy |ρ| < 1 for stationarity. + σ : float Standard deviation of the error term. y0 : float Initial value of the process at time t=0. @@ -152,11 +155,23 @@ class AR1(NamedTuple): T1 : int, optional Length of the future path to simulate (default is 100). """ - rho: float - sigma: float + ρ: float + σ: float y0: float - T0: int=100 - T1: int=100 + T0: int + T1: int + + +def make_ar1(ρ: float, σ: float, y0: float, T0: int = 100, T1: int = 100): + """ + Factory function to create an AR1 instance with default values for T0 and T1. + + Returns + ------- + AR1 + AR1 named tuple containing the specified parameters. + """ + return AR1(ρ=ρ, σ=σ, y0=y0, T0=T0, T1=T1) def AR1_simulate_past(ar1: AR1, key=key): @@ -166,7 +181,7 @@ def AR1_simulate_past(ar1: AR1, key=key): Parameters ---------- ar1 : AR1 - AR1 named tuple containing parameters (rho, sigma, y0, T0, T1). + AR1 named tuple containing parameters (ρ, σ, y0, T0, T1). key : jax.random.PRNGKey JAX random key for generating random noise. @@ -175,18 +190,18 @@ def AR1_simulate_past(ar1: AR1, key=key): initial_path : jax.numpy.ndarray Simulated path of the AR(1) process and the initial y0. """ - rho, sigma, y0, T0 = ar1.rho, ar1.sigma, ar1.y0, ar1.T0 - # Draw epsilons - eps = sigma * random.normal(key, (T0,)) + ρ, σ, y0, T0 = ar1.ρ, ar1.σ, ar1.y0, ar1.T0 + # Draw εs + ε = σ * random.normal(key, (T0,)) # Set step function - def ar1_step(y_prev, t_rho_eps): - rho, eps_t = t_rho_eps - y_t = rho * y_prev + eps_t + def ar1_step(y_prev, t_ρ_ε): + ρ, ε_t = t_ρ_ε + y_t = ρ * y_prev + ε_t return y_t, y_t # Scan over time steps - _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0, rho), eps)) + _, y_seq = lax.scan(ar1_step, y0, (jnp.full(T0, ρ), ε)) # Concatenate initial value initial_path = jnp.concatenate([jnp.array([y0]), y_seq]) @@ -204,7 +219,7 @@ def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): Parameters ---------- ar1 : AR1 - AR1 named tuple containing parameters (rho, sigma, y0, T0, T1). + AR1 named tuple containing parameters (ρ, σ, y0, T0, T1). y_T0 : float Value of the process at time T0. N: int @@ -217,16 +232,16 @@ def AR1_simulate_future(ar1: AR1, y_T0, N=10, key=key): future_path : jax.numpy.ndarray Simulated N paths of the AR(1) process of length T1. """ - rho, sigma, T1 = ar1.rho, ar1.sigma, ar1.T1 + ρ, σ, T1 = ar1.ρ, ar1.σ, ar1.T1 def single_path_scan(y_T0, subkey): - eps = sigma * random.normal(subkey, (T1,)) + ε = σ * random.normal(subkey, (T1,)) - def ar1_step(y_prev, t_rho_eps): - rho, eps_t = t_rho_eps - y_t = rho * y_prev + eps_t + def ar1_step(y_prev, t_ρ_ε): + ρ, ε_t = t_ρ_ε + y_t = ρ * y_prev + ε_t return y_t, y_t - _, y = lax.scan(ar1_step, y_T0, (jnp.full(T1, rho), eps)) + _, y = lax.scan(ar1_step, y_T0, (jnp.full(T1, ρ), ε)) return y # Split key to generate different paths @@ -249,7 +264,7 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): Parameters ---------- ar1 : AR1 - AR1 named tuple containing process parameters (rho, sigma, T0, T1). + AR1 named tuple containing process parameters (ρ, σ, T0, T1). initial_path : array-like Simulated initial path of the AR(1) process, shape(T0+1,). future_path : array-like @@ -266,13 +281,13 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): - 90% and 95% predictive confidence intervals - Expected future path """ - rho, sigma, T0, T1 = ar1.rho, ar1.sigma, ar1.T0, ar1.T1 + ρ, σ, T0, T1 = ar1.ρ, ar1.σ, ar1.T0, ar1.T1 # Compute moments and confidence intervals y_T0 = initial_path[-1] j = jnp.arange(1, T1+1) - center = rho**j * y_T0 - vars = sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) + center = ρ**j * y_T0 + vars = σ**2 * (1 - ρ**(2 * j)) / (1 - ρ**2) # 95% CI y_upper_c95 = center + 1.96 * jnp.sqrt(vars) @@ -313,11 +328,10 @@ def plot_path(ar1, initial_path, future_path, ax, key=key): --- mystnb: figure: - caption: | - Initial and predictive future paths + caption: "Initial and predictive future paths \n" name: fig_path --- -ar1 = AR1(rho=0.9, sigma=1, y0=10, T0=100, T1=100) +ar1 = make_ar1(ρ=0.9, σ=1, y0=10) # Simulate initial_path = AR1_simulate_past(ar1) @@ -438,24 +452,23 @@ Note that in defining the likelihood function, we choose to condition on the ini --- mystnb: figure: - caption: | - Posterior distributions (rho, sigma) - name: fig_post + caption: "AR(1) model" + name: fig_trace --- def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): """Draw a sample of size from the posterior distribution.""" def model(data): # Start with priors - rho = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable rho - sigma = numpyro.sample('sigma', dist.HalfNormal(jnp.sqrt(10))) + ρ = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable ρ + σ = numpyro.sample('sigma', dist.HalfNormal(jnp.sqrt(10))) # Define likelihood recursively for t in range(1, len(data)): # Expectation of y_t - mu = rho * data[t-1] + μ = ρ * data[t-1] # Likelihood of the actual realization. - numpyro.sample(f'y_{t}', dist.Normal(mu, sigma), obs=data[t]) + numpyro.sample(f'y_{t}', dist.Normal(μ, σ), obs=data[t]) # Compute posterior distribution of parameters nuts_kernel = NUTS(model) @@ -465,6 +478,7 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): nuts_kernel, num_warmup=5000, num_samples=size, + num_chains=4, # plot 4 chains in the trace progress_bar=False) # Run MCMC @@ -476,24 +490,26 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): 'sigma': mcmc.get_samples()['sigma'] } - # Plot posterior distributions + # Plot posterior distributions and trace plots if dis_plot == 1: - fig, ax = plt.subplots(1, 2, figsize=(12, 5)) - sns.histplot( - post_sample['rho'], kde=True, stat="density", bins=bins, ax=ax[0] + plot_data = az.from_numpyro(posterior=mcmc) + axes = az.plot_trace( + data=plot_data, + compact=True, + lines=[ + ("ρ", {}, ar1.ρ), + ("σ", {}, ar1.σ), + ], + backend_kwargs={"figsize": (10, 6), "layout": "constrained"}, ) - ax[0].set_xlabel(r"$\rho$") - sns.histplot( - post_sample['sigma'], kde=True, stat="density", bins=bins, ax=ax[1] - ) - ax[1].set_xlabel(r"$\sigma$") + return post_sample post_samples = draw_from_posterior(initial_path) ``` -The graphs above portray posterior distributions. +The graphs above portray posterior distributions and trace plots. The posterior distributions (top row) show the marginal distributions of the parameters after observing the data, while the trace plots (bottom row) help diagnose MCMC convergence by showing how the sampler explored the parameter space over iterations. ## Calculating Sample Path Statistics @@ -567,6 +583,7 @@ def compute_path_statistics(initial_path, future_path): ) return path_stats ``` + The following function creates visualizations of the path statistics in a subplot grid. ```{code-cell} ipython3 @@ -611,7 +628,7 @@ def plot_Wecker(ar1: AR1, initial_path, ax, N=1000): Parameters ---------- ar1 : AR1 - An AR1 named tuple containing the process parameters (rho, sigma, T0, T1). + An AR1 named tuple containing the process parameters (ρ, σ, T0, T1). initial_path : array-like The initial observed path of the AR(1) process. N : int @@ -669,14 +686,14 @@ def plot_extended_Wecker( index = random.choice( key, jnp.arange(len(post_samples['rho'])), (N + 1,), replace=False ) - rho_sample = post_samples['rho'][index] - sigma_sample = post_samples['sigma'][index] + ρ_sample = post_samples['rho'][index] + σ_sample = post_samples['sigma'][index] # Compute path statistics subkeys = random.split(key, num=N) def step(carry, n): - ar1_n = AR1(rho=rho_sample[n], sigma=sigma_sample[n], y0=y0) + ar1_n = make_ar1(ρ=ρ_sample[n], σ=σ_sample[n], y0=y0, T1=T1) future_temp = AR1_simulate_future( ar1_n, y_T0, N=1, key=subkeys[n] ).reshape(-1) From 56d96a968f9fb1c77000a3433bd4968a5305b3bf Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Fri, 5 Sep 2025 16:41:37 +0800 Subject: [PATCH 09/16] set numpyro.set_host_device_count --- lectures/ar1_turningpts.md | 1 + 1 file changed, 1 insertion(+) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index af65a0862..16ae684c8 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -66,6 +66,7 @@ import arviz as az sns.set_style('white') colors = sns.color_palette() key = random.PRNGKey(0) +numpyro.set_host_device_count(4) ``` ## A Univariate First-Order Autoregressive Process From d668cf4edacafc7ec0576939f061d1a0f4dfd82e Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Fri, 5 Sep 2025 17:45:02 +0800 Subject: [PATCH 10/16] modify var_nam in arviz --- lectures/ar1_turningpts.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 16ae684c8..102310afe 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -498,8 +498,8 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): data=plot_data, compact=True, lines=[ - ("ρ", {}, ar1.ρ), - ("σ", {}, ar1.σ), + ("rho", {}, ar1.ρ), + ("sigma", {}, ar1.σ), ], backend_kwargs={"figsize": (10, 6), "layout": "constrained"}, ) From 33735e2dabcc2d96dc8252c784a061908f1b78c9 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Sun, 7 Sep 2025 13:23:25 +0800 Subject: [PATCH 11/16] Reorder numpyro imports and device count setup --- lectures/ar1_turningpts.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 102310afe..c3d2316f1 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -49,24 +49,24 @@ import matplotlib.pyplot as plt import seaborn as sns from typing import NamedTuple +# numpyro +import numpyro +import numpyro.distributions as dist +from numpyro.infer import MCMC, NUTS +numpyro.set_host_device_count(4) + # jax import jax import jax.random as random import jax.numpy as jnp from jax import lax -# numpyro -import numpyro -import numpyro.distributions as dist -from numpyro.infer import MCMC, NUTS - # arviz import arviz as az sns.set_style('white') colors = sns.color_palette() key = random.PRNGKey(0) -numpyro.set_host_device_count(4) ``` ## A Univariate First-Order Autoregressive Process From 2a37d97c743943d36c89b4dbda0016e2e2bbcade Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Sun, 7 Sep 2025 14:19:18 +0800 Subject: [PATCH 12/16] test local_device_count --- lectures/ar1_turningpts.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index c3d2316f1..c2b78bb29 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -67,6 +67,9 @@ import arviz as az sns.set_style('white') colors = sns.color_palette() key = random.PRNGKey(0) + +print(f"Available devices: {jax.local_device_count()}") +print(f"Device list: {jax.devices()}") ``` ## A Univariate First-Order Autoregressive Process From 986b44578817e73ee5b0c401df93744103b6598d Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 8 Sep 2025 15:53:16 +1000 Subject: [PATCH 13/16] update the vectorization strategy --- lectures/ar1_turningpts.md | 48 ++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index c2b78bb29..a04e4b3b2 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.17.1 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -53,7 +53,6 @@ from typing import NamedTuple import numpyro import numpyro.distributions as dist from numpyro.infer import MCMC, NUTS -numpyro.set_host_device_count(4) # jax import jax @@ -67,9 +66,6 @@ import arviz as az sns.set_style('white') colors = sns.color_palette() key = random.PRNGKey(0) - -print(f"Available devices: {jax.local_device_count()}") -print(f"Device list: {jax.devices()}") ``` ## A Univariate First-Order Autoregressive Process @@ -456,16 +452,17 @@ Note that in defining the likelihood function, we choose to condition on the ini --- mystnb: figure: - caption: "AR(1) model" + caption: AR(1) model name: fig_trace --- def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): - """Draw a sample of size from the posterior distribution.""" + """Draw a sample of size from the posterior distribution.""" + def model(data): # Start with priors - ρ = numpyro.sample('rho', dist.Uniform(-1, 1)) # Assume stable ρ - σ = numpyro.sample('sigma', dist.HalfNormal(jnp.sqrt(10))) - + ρ = numpyro.sample('ρ', dist.Uniform(-1, 1)) # Assume stable ρ + σ = numpyro.sample('σ', dist.HalfNormal(jnp.sqrt(10))) + # Define likelihood recursively for t in range(1, len(data)): # Expectation of y_t @@ -476,24 +473,26 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): # Compute posterior distribution of parameters nuts_kernel = NUTS(model) - + # Define MCMC class to compute the posteriors mcmc = MCMC( nuts_kernel, num_warmup=5000, num_samples=size, - num_chains=4, # plot 4 chains in the trace - progress_bar=False) - + num_chains=4, # plot 4 chains in the trace + progress_bar=False, + chain_method='vectorized' + ) + # Run MCMC mcmc.run(key, data=data) - + # Get posterior samples post_sample = { - 'rho': mcmc.get_samples()['rho'], - 'sigma': mcmc.get_samples()['sigma'] + 'ρ': mcmc.get_samples()['ρ'], + 'σ': mcmc.get_samples()['σ'], } - + # Plot posterior distributions and trace plots if dis_plot == 1: plot_data = az.from_numpyro(posterior=mcmc) @@ -501,14 +500,13 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): data=plot_data, compact=True, lines=[ - ("rho", {}, ar1.ρ), - ("sigma", {}, ar1.σ), + ("σ", {}, ar1.ρ), + ("ρ", {}, ar1.σ), ], backend_kwargs={"figsize": (10, 6), "layout": "constrained"}, ) - - return post_sample + return post_sample post_samples = draw_from_posterior(initial_path) ``` @@ -688,10 +686,10 @@ def plot_extended_Wecker( # Select a parameter sample index = random.choice( - key, jnp.arange(len(post_samples['rho'])), (N + 1,), replace=False + key, jnp.arange(len(post_samples['ρ'])), (N + 1,), replace=False ) - ρ_sample = post_samples['rho'][index] - σ_sample = post_samples['sigma'][index] + ρ_sample = post_samples['ρ'][index] + σ_sample = post_samples['σ'][index] # Compute path statistics subkeys = random.split(key, num=N) From 6de7ab1e34e42cd18e9ebd9cfa5fcb440eab515b Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 8 Sep 2025 14:26:02 +0800 Subject: [PATCH 14/16] fix typo and add explanation --- lectures/ar1_turningpts.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index a04e4b3b2..8246846de 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -172,8 +172,11 @@ def make_ar1(ρ: float, σ: float, y0: float, T0: int = 100, T1: int = 100): AR1 named tuple containing the specified parameters. """ return AR1(ρ=ρ, σ=σ, y0=y0, T0=T0, T1=T1) +``` +Using the `AR1` class, we can simulate paths more conveniently. The following function simulates an initial path with $T0$ length. +```{code-cell} ipython3 def AR1_simulate_past(ar1: AR1, key=key): """ Simulate a realization of the AR(1) process for T0 periods. @@ -500,8 +503,8 @@ def draw_from_posterior(data, size=10000, bins=20, dis_plot=1, key=key): data=plot_data, compact=True, lines=[ - ("σ", {}, ar1.ρ), - ("ρ", {}, ar1.σ), + ("ρ", {}, ar1.ρ), + ("σ", {}, ar1.σ), ], backend_kwargs={"figsize": (10, 6), "layout": "constrained"}, ) From 0376354d9fae89c6d0481c76e85f080e67319347 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 23 Sep 2025 15:22:24 +0900 Subject: [PATCH 15/16] address review comments --- lectures/ar1_turningpts.md | 64 ++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 31 deletions(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index 8246846de..aa99d3c4e 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -158,20 +158,8 @@ class AR1(NamedTuple): ρ: float σ: float y0: float - T0: int - T1: int - - -def make_ar1(ρ: float, σ: float, y0: float, T0: int = 100, T1: int = 100): - """ - Factory function to create an AR1 instance with default values for T0 and T1. - - Returns - ------- - AR1 - AR1 named tuple containing the specified parameters. - """ - return AR1(ρ=ρ, σ=σ, y0=y0, T0=T0, T1=T1) + T0: int = 100 + T1: int = 100 ``` Using the `AR1` class, we can simulate paths more conveniently. The following function simulates an initial path with $T0$ length. @@ -334,7 +322,7 @@ mystnb: caption: "Initial and predictive future paths \n" name: fig_path --- -ar1 = make_ar1(ρ=0.9, σ=1, y0=10) +ar1 = AR1(ρ=0.9, σ=1, y0=10) # Simulate initial_path = AR1_simulate_past(ar1) @@ -346,9 +334,7 @@ plot_path(ar1, initial_path, future_path, ax) plt.show() ``` -As functions of forecast horizon, the coverage intervals have shapes like those described in -https://python.quantecon.org/perm_income_cons.html - +As functions of forecast horizon, the coverage intervals have shapes like those described in [Permanent Income II: LQ Techniques](perm_income_cons) ## Predictive Distributions of Path Properties @@ -644,16 +630,24 @@ def plot_Wecker(ar1: AR1, initial_path, ax, N=1000): future_path = AR1_simulate_future(ar1, y_T0, N=N) plot_path(ar1, initial_path, future_path, ax[0, 0]) + next_reces = jnp.zeros(N) + severe_rec = jnp.zeros(N) + min_val_8q = jnp.zeros(N) + next_up_turn = jnp.zeros(N) + next_down_turn = jnp.zeros(N) + # Simulate future paths and compute statistics - def step(carry, n): + for n in range(N): future_temp = future_path[n, :] - (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn + (next_reces_val, severe_rec_val, min_val_8q_val, + next_up_turn_val, next_down_turn_val ) = compute_path_statistics(initial_path, future_temp) - return carry, (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn) - - _, (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn - ) = lax.scan(step, None, jnp.arange(N)) + next_reces = next_reces.at[n].set(next_reces_val) + severe_rec = severe_rec.at[n].set(severe_rec_val) + min_val_8q = min_val_8q.at[n].set(min_val_8q_val) + next_up_turn = next_up_turn.at[n].set(next_up_turn_val) + next_down_turn = next_down_turn.at[n].set(next_down_turn_val) # Plot path statistics plot_path_stats(next_reces, severe_rec, min_val_8q, @@ -695,20 +689,28 @@ def plot_extended_Wecker( σ_sample = post_samples['σ'][index] # Compute path statistics + next_reces = jnp.zeros(N) + severe_rec = jnp.zeros(N) + min_val_8q = jnp.zeros(N) + next_up_turn = jnp.zeros(N) + next_down_turn = jnp.zeros(N) + subkeys = random.split(key, num=N) - def step(carry, n): - ar1_n = make_ar1(ρ=ρ_sample[n], σ=σ_sample[n], y0=y0, T1=T1) + for n in range(N): + ar1_n = AR1(ρ=ρ_sample[n], σ=σ_sample[n], y0=y0, T1=T1) future_temp = AR1_simulate_future( ar1_n, y_T0, N=1, key=subkeys[n] ).reshape(-1) - (next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn + (next_reces_val, severe_rec_val, min_val_8q_val, + next_up_turn_val, next_down_turn_val ) = compute_path_statistics(initial_path, future_temp) - return carry, (future_temp, next_reces, severe_rec, - min_val_8q, next_up_turn, next_down_turn) - _, (future_path, next_reces, severe_rec, min_val_8q, next_up_turn, next_down_turn - ) = jax.lax.scan(step, None, jnp.arange(N)) + next_reces = next_reces.at[n].set(next_reces_val) + severe_rec = severe_rec.at[n].set(severe_rec_val) + min_val_8q = min_val_8q.at[n].set(min_val_8q_val) + next_up_turn = next_up_turn.at[n].set(next_up_turn_val) + next_down_turn = next_down_turn.at[n].set(next_down_turn_val) # Plot simulated initial and future paths plot_path(ar1, initial_path, future_path, ax[0, 0]) From f802f5165d5e18bd6ab3b05a7199379fabfff7f3 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 23 Sep 2025 16:23:25 +0900 Subject: [PATCH 16/16] fix document refrence --- lectures/ar1_turningpts.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/ar1_turningpts.md b/lectures/ar1_turningpts.md index aa99d3c4e..e5360cdc2 100644 --- a/lectures/ar1_turningpts.md +++ b/lectures/ar1_turningpts.md @@ -334,7 +334,7 @@ plot_path(ar1, initial_path, future_path, ax) plt.show() ``` -As functions of forecast horizon, the coverage intervals have shapes like those described in [Permanent Income II: LQ Techniques](perm_income_cons) +As functions of forecast horizon, the coverage intervals have shapes like those described in {doc}`perm_income_cons` ## Predictive Distributions of Path Properties