- prospective values $y_{t+j}$ of a random process $\{y_t\}$ that is governed by the AR(1) process

- sample path properties that are defined as non-linear functions of future values $\{y_{t+j}\}_{j \geq 1}$ at time $t$.

- sample path properties that are defined as non-linear functions of future values $\{y_{t+j}\}_{j \geq 1}$ at time $t$

**Sample path properties** are things like "time to next turning point" or "time to next recession"

**Sample path properties** are things like "time to next turning point" or "time to next recession".

To investigate sample path properties we'll use a simulation procedure recommended by Wecker {cite}`wecker1979predicting`.

```{code-cell} ipython3

import numpy as np

import arviz as az

import pymc as pmc

import matplotlib.pyplot as plt

import seaborn as sns

sns.set_style('white')

colors = sns.color_palette()

import logging

logging.basicConfig()

logger = logging.getLogger('pymc')

- the time until the next "recession",

- the minimum value of $Y$ over the next 8 periods,

- "severe recession", and

- the time until the next turning point (positive or negative)

- the time until the next turning point (positive or negative).

To accomplish that for situations in which we are uncertain about parameter values, we shall extend the approach Wecker {cite}`wecker1979predicting` in the following way.

To accomplish that for situations in which we are uncertain about parameter values, we shall extend Wecker's {cite}`wecker1979predicting` approach in the following way.

- first simulate an initial path of length $T_0$;

- for a given prior, draw a sample of size $N$ from the posterior joint distribution of parameters $\left(\rho,\sigma\right)$ after observing the initial path;

# Allocate space and draw epsilons

y = np.empty(T)

eps = np.random.normal(0,sigma,T)

eps = np.random.normal(0, sigma, T)

# Initial condition and step forward

y[0] = y0

for t in range(1, T):

y[t] = rho*y[t-1] + eps[t]

y[t] = rho * y[t-1] + eps[t]

return y

"""

Plot the initial path and the preceding predictive densities

"""

# compute .9 confidence interval]

# Compute .9 confidence interval]

y0 = initial_path[-1]

center = np.array([rho**j*y0 for j in range(T1)])

vars = np.array([sigma**2*(1-rho**(2*j))/(1-rho**2) for j in range(T1)])

center = np.array([rho**j * y0 for j in range(T1)])

vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])

y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)

y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)

# plot

fig, ax = plt.subplots(1,1, figsize = (12, 6))

# Plot

fig, ax = plt.subplots(1, 1, figsize=(12, 6))

ax.set_title("Initial Path and Predictive Densities", fontsize=15)

ax.plot(np.arange(-T0+1, 1), initial_path)

ax.plot(np.arange(-T0 + 1, 1), initial_path)

ax.set_xlim([-T0, T1])

ax.axvline(0, linestyle='--', alpha=.4, color='k', lw=1)

# simulate future paths

# Simulate future paths

for i in range(10):

y_future = AR1_simulate(rho, sigma, y0, T1)

ax.plot(np.arange(T1), y_future, color='grey', alpha=.5)

# plot 90% CI

# Plot 90% CI

ax.fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3, label='95% CI')

ax.fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35, label='90% CI')

ax.plot(np.arange(T1), center, color='red', alpha=.7, label='expected mean')

T0, T1 = 100, 100

y0 = 10

# simulate

# Simulate

np.random.seed(145)

initial_path = AR1_simulate(rho, sigma, y0, T0)

# plot

# Plot

plot_initial_path(initial_path)

```

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$.

$$

Wecker {cite}`wecker1979predicting` also studied **the minimum value of $Y$ over the next 8 quarters**

which can be defined as the random variable

which can be defined as the random variable.

$$

M_t(\omega) := \min \{ Y_{t+1}(\omega); Y_{t+2}(\omega); \dots; Y_{t+8}(\omega)\}

with AR1_model:

#start with priors

rho = pmc.Uniform('rho',lower=-1.,upper=1.) #assume stable rho

# Start with priors

rho = pmc.Uniform('rho',lower=-1.,upper=1.) # Assume stable rho

sigma = pmc.HalfNormal('sigma', sigma = np.sqrt(10))

# Expected value of y at the next period (rho * y)

yhat = rho*sample[:-1]

yhat = rho * sample[:-1]

# Likelihood of the actual realization.

y_like = pmc.Normal('y_obs', mu = yhat, sigma=sigma, observed=sample[1:])

y_like = pmc.Normal('y_obs', mu=yhat, sigma=sigma, observed=sample[1:])

with AR1_model:

trace = pmc.sample(10000, tune=5000)

# check condition

with AR1_model:

az.plot_trace(trace, figsize=(17,6))

az.plot_trace(trace, figsize=(17, 6))

rhos = trace.posterior.rho.values.flatten()

sigmas = trace.posterior.sigma.values.flatten()

## Calculating Sample Path Statistics

Our next step is to prepare Python codeto compute our sample path statistics.

Our next step is to prepare Python code to compute our sample path statistics.

```{code-cell} ipython3

# define statistics

n = z.shape[0]

sr = (z < -.02).astype(int)

indices = np.where(sr==1)[0]

indices = np.where(sr == 1)[0]

if len(indices)==0:

if len(indices) == 0:

return T1

else:

return indices[0] + 1

(omega[i+2] > omega[i+3]) and (omega[i+3] > omega[i+4])):

T[i] = -1

up_turn = np.where(T==1)[0][0] + 1 if (1 in T) == True else T1

down_turn = np.where(T==-1)[0][0] + 1 if (-1 in T) == True else T1

up_turn = np.where(T == 1)[0][0] + 1 if (1 in T) == True else T1

down_turn = np.where(T == -1)[0][0] + 1 if (-1 in T) == True else T1

return up_turn, down_turn

```

"""

Plot the predictive distributions from "pure" Wecker's method.

"""

# store outcomes

# Store outcomes

next_reces = np.zeros(N)

severe_rec = np.zeros(N)

min_vals = np.zeros(N)

next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)

# compute .9 confidence interval]

# Compute .9 confidence interval]

y0 = initial_path[-1]

center = np.array([rho**j*y0 for j in range(T1)])

vars = np.array([sigma**2*(1-rho**(2*j))/(1-rho**2) for j in range(T1)])

center = np.array([rho**j * y0 for j in range(T1)])

vars = np.array([sigma**2 * (1 - rho**(2 * j)) / (1 - rho**2) for j in range(T1)])

y_bounds1_c95, y_bounds2_c95 = center + 1.96 * np.sqrt(vars), center - 1.96 * np.sqrt(vars)

y_bounds1_c90, y_bounds2_c90 = center + 1.65 * np.sqrt(vars), center - 1.65 * np.sqrt(vars)

# plot

# Plot

ax[0, 0].set_title("Initial path and predictive densities", fontsize=15)

ax[0, 0].plot(np.arange(-T0+1, 1), initial_path)

ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)

ax[0, 0].set_xlim([-T0, T1])

ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)

# plot 90% CI

# Plot 90% CI

ax[0, 0].fill_between(np.arange(T1), y_bounds1_c95, y_bounds2_c95, alpha=.3)

ax[0, 0].fill_between(np.arange(T1), y_bounds1_c90, y_bounds2_c90, alpha=.35)

ax[0, 0].plot(np.arange(T1), center, color='red', alpha=.7)

# simulate future paths

# Simulate future paths

for n in range(N):

sim_path = AR1_simulate(rho, sigma, initial_path[-1], T1)

next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))

if n%(N/10) == 0:

ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)

# return next_up_turn, next_down_turn

# Return next_up_turn, next_down_turn

sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.8, label='True parameters')

ax[0, 1].set_title("Predictive distribution of time until the next recession", fontsize=13)

Now we apply we apply our "extended" Wecker method based on predictive densities of $y$ defined by

{eq}`ar1-tp-eq4` that acknowledge posterior uncertainty in the parameters $\rho, \sigma$.

To approximate the intergration on the right side of {eq}`ar1-tp-eq4`, we repeately draw parameters from the joint posterior distribution each time we simulate a sequence of future values from model {eq}`ar1-tp-eq1`.

To approximate the intergration on the right side of {eq}`ar1-tp-eq4`, we repeatedly draw parameters from the joint posterior distribution each time we simulate a sequence of future values from model {eq}`ar1-tp-eq1`.

```{code-cell} ipython3

def plot_extended_Wecker(post_samples, initial_path, N, ax):

"""

Plot the extended Wecker's predictive distribution

"""

# select a sample

index = np.random.choice(np.arange(len(post_samples['rho'])), N+1, replace=False)

# Select a sample

index = np.random.choice(np.arange(len(post_samples['rho'])), N + 1, replace=False)

rho_sample = post_samples['rho'][index]

sigma_sample = post_samples['sigma'][index]

# store outcomes

# Store outcomes

next_reces = np.zeros(N)

severe_rec = np.zeros(N)

min_vals = np.zeros(N)

next_up_turn, next_down_turn = np.zeros(N), np.zeros(N)

# plot

# Plot

ax[0, 0].set_title("Initial path and future paths simulated from posterior draws", fontsize=15)

ax[0, 0].plot(np.arange(-T0+1, 1), initial_path)

ax[0, 0].plot(np.arange(-T0 + 1, 1), initial_path)

ax[0, 0].set_xlim([-T0, T1])

ax[0, 0].axvline(0, linestyle='--', alpha=.4, color='k', lw=1)

# simulate future paths

# Simulate future paths

for n in range(N):

sim_path = AR1_simulate(rho_sample[n], sigma_sample[n], initial_path[-1], T1)

next_reces[n] = next_recession(np.hstack([initial_path[-3:-1], sim_path]))

severe_rec[n] = severe_recession(sim_path)

min_vals[n] = minimum_value(sim_path)

next_up_turn[n], next_down_turn[n] = next_turning_point(sim_path)

if n%(N/10) == 0:

if n % (N / 10) == 0:

ax[0, 0].plot(np.arange(T1), sim_path, color='gray', alpha=.3, lw=1)

# return next_up_turn, next_down_turn

# Return next_up_turn, next_down_turn

sns.histplot(next_reces, kde=True, stat='density', ax=ax[0, 1], alpha=.6, color=colors[1], label='Sampling from posterior')

ax[0, 1].set_title("Predictive distribution of time until the next recession", fontsize=13)

sns.histplot(next_down_turn, kde=True, stat='density', ax=ax[2, 1], alpha=.6, color=colors[1], label='Sampling from posterior')

ax[2, 1].set_title("Predictive distribution of time until the next negative turn", fontsize=13)

fig, ax = plt.subplots(3, 2, figsize=(15,12))

fig, ax = plt.subplots(3, 2, figsize=(15, 12))

plot_extended_Wecker(post_samples, initial_path, 1000, ax)

plt.show()

```

## Comparison

Finally, we plot both the original Wecker method and the extended method with parameter values drawn from the posterior together to compare the differnces that emerge from pretending to know parameter values when they are actually uncertain.

Finally, we plot both the original Wecker method and the extended method with parameter values drawn from the posterior together to compare the differences that emerge from pretending to know parameter values when they are actually uncertain.

```{code-cell} ipython3

fig, ax = plt.subplots(3, 2, figsize=(15,12))

plot_Wecker(initial_path, 1000, ax)

ax[0,0].clear()

ax[0, 0].clear()

plot_extended_Wecker(post_samples, initial_path, 1000, ax)

plt.legend()

plt.show()