diff --git a/lectures/divergence_measures.md b/lectures/divergence_measures.md index b95d30d4e..747b815b3 100644 --- a/lectures/divergence_measures.md +++ b/lectures/divergence_measures.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.1 + jupytext_version: 1.16.6 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -60,10 +60,6 @@ import pandas as pd from IPython.display import display, Math ``` - - - - ## Primer on entropy, cross-entropy, KL divergence Before diving in, we'll introduce some useful concepts in a simple setting. @@ -163,6 +159,8 @@ f(z; a, b) = \frac{\Gamma(a+b) z^{a-1} (1-z)^{b-1}}{\Gamma(a) \Gamma(b)} \Gamma(p) := \int_{0}^{\infty} x^{p-1} e^{-x} dx $$ +We introduce two Beta distributions $f(x)$ and $g(x)$, which we will use to illustrate the different divergence measures. + Let's define parameters and density functions in Python ```{code-cell} ipython3 @@ -198,8 +196,6 @@ plt.legend() plt.show() ``` - - (rel_entropy)= ## Kullback–Leibler divergence @@ -457,15 +453,14 @@ plt.show() We now generate plots illustrating how overlap visually diminishes as divergence measures increase. - ```{code-cell} ipython3 param_grid = [ ((1, 1), (1, 1)), ((1, 1), (1.5, 1.2)), ((1, 1), (2, 1.5)), ((1, 1), (3, 1.2)), - ((1, 1), (5, 1)), - ((1, 1), (0.3, 0.3)) + ((1, 1), (0.3, 0.3)), + ((1, 1), (5, 1)) ] ``` @@ -512,8 +507,6 @@ def plot_dist_diff(para_grid): divergence_data = plot_dist_diff(param_grid) ``` - - ## KL divergence and maximum-likelihood estimation diff --git a/lectures/imp_sample.md b/lectures/imp_sample.md index 11c7258f9..09ae67576 100644 --- a/lectures/imp_sample.md +++ b/lectures/imp_sample.md @@ -36,7 +36,7 @@ import matplotlib.pyplot as plt from math import gamma ``` -## Mathematical Expectation of Likelihood Ratio +## Mathematical expectation of likelihood ratio In {doc}`this lecture `, we studied a likelihood ratio $\ell \left(\omega_t\right)$ @@ -57,11 +57,10 @@ $$ Our goal is to approximate the mathematical expectation $E \left[ L\left(\omega^t\right) \right]$ well. In {doc}`this lecture `, we showed that $E \left[ L\left(\omega^t\right) \right]$ equals $1$ for all $t$. + We want to check out how well this holds if we replace $E$ by with sample averages from simulations. -This turns out to be easier said than done because for -Beta distributions assumed above, $L\left(\omega^t\right)$ has -a very skewed distribution with a very long tail as $t \rightarrow \infty$. +This turns out to be easier said than done because for Beta distributions assumed above, $L\left(\omega^t\right)$ has a very skewed distribution with a very long tail as $t \rightarrow \infty$. This property makes it difficult efficiently and accurately to estimate the mean by standard Monte Carlo simulation methods. @@ -156,7 +155,7 @@ $$ E^g\left[\ell\left(\omega\right)\right] = \int_\Omega \ell(\omega) g(\omega) d\omega = \int_\Omega \ell(\omega) \frac{g(\omega)}{h(\omega)} h(\omega) d\omega = E^h\left[\ell\left(\omega\right) \frac{g(\omega)}{h(\omega)}\right] $$ -## Selecting a Sampling Distribution +## Selecting a sampling distribution Since we must use an $h$ that has larger mass in parts of the distribution to which $g$ puts low mass, we use $h=Beta(0.5, 0.5)$ as our importance distribution. @@ -178,7 +177,7 @@ plt.ylim([0., 3.]) plt.show() ``` -## Approximating a Cumulative Likelihood Ratio +## Approximating a cumulative likelihood ratio We now study how to use importance sampling to approximate ${E} \left[L(\omega^t)\right] = \left[\prod_{i=1}^T \ell \left(\omega_i\right)\right]$. @@ -233,7 +232,7 @@ For our importance sampling estimate, we set $q = h$. estimate(g_a, g_b, h_a, h_b, T=1, N=10000) ``` -Evidently, even at T=1, our importance sampling estimate is closer to $1$ than is the Monte Carlo estimate. +Evidently, even at $T=1$, our importance sampling estimate is closer to $1$ than is the Monte Carlo estimate. Bigger differences arise when computing expectations over longer sequences, $E_0\left[L\left(\omega^t\right)\right]$. @@ -248,7 +247,17 @@ estimate(g_a, g_b, g_a, g_b, T=10, N=10000) estimate(g_a, g_b, h_a, h_b, T=10, N=10000) ``` -## Distribution of Sample Mean +The Monte Carlo method underestimates because the likelihood ratio $L(\omega^T) = \prod_{t=1}^T \frac{f(\omega_t)}{g(\omega_t)}$ has a highly skewed distribution under $g$. + +Most samples from $g$ produce small likelihood ratios, while the true mean requires occasional very large values that are rarely sampled. + +In our case, since $g(\omega) \to 0$ as $\omega \to 0$ while $f(\omega)$ remains constant, the Monte Carlo procedure undersamples precisely where the likelihood ratio $\frac{f(\omega)}{g(\omega)}$ is largest. + +As $T$ increases, this problem worsens exponentially, making standard Monte Carlo increasingly unreliable. + +Importance sampling with $q = h$ fixes this by sampling more uniformly from regions important to both $f$ and $g$. + +## Distribution of sample mean We next study the bias and efficiency of the Monte Carlo and importance sampling approaches. @@ -323,7 +332,7 @@ The simulation exercises above show that the importance sampling estimates are u Evidently, the bias increases with increases in $T$. -## Choosing a Sampling Distribution +## Choosing a sampling distribution +++ @@ -331,9 +340,7 @@ Above, we arbitraily chose $h = Beta(0.5,0.5)$ as the importance distribution. Is there an optimal importance distribution? -In our particular case, since we know in advance that $E_0 \left[ L\left(\omega^t\right) \right] = 1$. - -We can use that knowledge to our advantage. +In our particular case, since we know in advance that $E_0 \left[ L\left(\omega^t\right) \right] = 1$, we can use that knowledge to our advantage. Thus, suppose that we simply use $h = f$. @@ -364,10 +371,10 @@ b_list = [0.5, 1.2, 5.] ```{code-cell} ipython3 w_range = np.linspace(1e-5, 1-1e-5, 1000) -plt.plot(w_range, g(w_range), label=f'p=Beta({g_a}, {g_b})') -plt.plot(w_range, p(w_range, a_list[0], b_list[0]), label=f'g=Beta({a_list[0]}, {b_list[0]})') -plt.plot(w_range, p(w_range, a_list[1], b_list[1]), label=f'g=Beta({a_list[1]}, {b_list[1]})') -plt.plot(w_range, p(w_range, a_list[2], b_list[2]), label=f'g=Beta({a_list[2]}, {b_list[2]})') +plt.plot(w_range, g(w_range), label=f'g=Beta({g_a}, {g_b})') +plt.plot(w_range, p(w_range, a_list[0], b_list[0]), label=f'$h_1$=Beta({a_list[0]},{b_list[0]})') +plt.plot(w_range, p(w_range, a_list[1], b_list[1]), label=f'$h_2$=Beta({a_list[1]},{b_list[1]})') +plt.plot(w_range, p(w_range, a_list[2], b_list[2]), label=f'$h_3$=Beta({a_list[2]},{b_list[2]})') plt.title('real data generating process $g$ and importance distribution $h$') plt.legend() plt.ylim([0., 3.]) diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index 5a7213702..fc62a4145 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -26,6 +26,12 @@ kernelspec: :depth: 2 ``` +In addition to what’s in Anaconda, this lecture will need the following libraries: + +```{code-cell} ipython3 +:tags: [hide-output] +!pip install --upgrade quantecon +``` ## Overview @@ -194,7 +200,7 @@ def plot_likelihood_paths(l_seq, title="Likelihood ratio paths", ``` (nature_likeli)= -## Nature Permanently Draws from Density g +## Nature permanently draws from density g We first simulate the likelihood ratio process when nature permanently draws from $g$. @@ -264,7 +270,7 @@ Mathematical induction implies $E\left[L\left(w^{t}\right)\bigm|q=g\right]=1$ for all $t \geq 1$. -## Peculiar Property +## Peculiar property How can $E\left[L\left(w^{t}\right)\bigm|q=g\right]=1$ possibly be true when most probability mass of the likelihood ratio process is piling up near $0$ as @@ -300,7 +306,7 @@ We explain the problem in more detail in {doc}`this lecture `. There we describe an alternative way to compute the mean of a likelihood ratio by computing the mean of a _different_ random variable by sampling from a _different_ probability distribution. -## Nature Permanently Draws from Density f +## Nature permanently draws from density f Now suppose that before time $0$ nature permanently decided to draw repeatedly from density $f$. @@ -348,7 +354,7 @@ plt.plot(range(T), np.sum(l_seq_f > 10000, axis=0) / N) plt.show() ``` -## Likelihood Ratio Test +## Likelihood ratio test We now describe how to employ the machinery of Neyman and Pearson {cite}`Neyman_Pearson` to test the hypothesis that history $w^t$ is generated by repeated @@ -959,7 +965,7 @@ $$ p(\textrm{wrong decision}) = {1 \over 2} (\alpha_T + \beta_T) . $$ (eq:detectionerrorprob) -Now let's simulate timing protocol 1 and 2 and compute the error probabilities +Now let's simulate timing protocol 1 and compute the error probabilities ```{code-cell} ipython3 @@ -993,33 +999,9 @@ def compute_protocol_1_errors(π_minus_1, T_max, N_simulations, f_func, g_func, 'L_cumulative': L_cumulative, 'true_models': true_models } - -def compute_protocol_2_errors(π_minus_1, T_max, N_simulations, f_func, g_func, - F_params=(1, 1), G_params=(3, 1.2)): - """ - Compute error probabilities for Protocol 2. - """ - sequences, true_models = protocol_2(π_minus_1, - T_max, N_simulations, F_params, G_params) - l_ratios, _ = compute_likelihood_ratios(sequences, f_func, g_func) - - T_range = np.arange(1, T_max + 1) - - accuracy = np.empty(T_max) - for t in range(T_max): - predictions = (l_ratios[:, t] >= 1) - actual = true_models[:, t] - accuracy[t] = np.mean(predictions == actual) - - return { - 'T_range': T_range, - 'accuracy': accuracy, - 'l_ratios': l_ratios, - 'true_models': true_models - } ``` -The following code visualizes the error probabilities for timing protocol 1 and 2 +The following code visualizes the error probabilities for timing protocol 1 ```{code-cell} ipython3 :tags: [hide-input] @@ -1058,44 +1040,6 @@ def analyze_protocol_1(π_minus_1, T_max, N_simulations, f_func, g_func, return result -def analyze_protocol_2(π_minus_1, T_max, N_simulations, f_func, g_func, - theory_error=None, F_params=(1, 1), G_params=(3, 1.2)): - """Analyze Protocol 2.""" - result = compute_protocol_2_errors(π_minus_1, T_max, N_simulations, - f_func, g_func, F_params, G_params) - - # Plot results - plt.figure(figsize=(10, 6)) - plt.plot(result['T_range'], result['accuracy'], - 'b-', linewidth=2, label='empirical accuracy') - - if theory_error is not None: - plt.axhline(1 - theory_error, color='r', linestyle='--', - label=f'theoretical accuracy = {1 - theory_error:.4f}') - - plt.xlabel('$t$') - plt.ylabel('accuracy') - plt.legend() - plt.ylim(0.5, 1.0) - plt.show() - - return result - -def compare_protocols(result1, result2): - """Compare results from both protocols.""" - plt.figure(figsize=(10, 6)) - - plt.plot(result1['T_range'], result1['error_prob'], linewidth=2, - label='Protocol 1 (Model Selection)') - plt.plot(result2['T_range'], 1 - result2['accuracy'], - linestyle='--', linewidth=2, - label='Protocol 2 (classification)') - - plt.xlabel('$T$') - plt.ylabel('error probability') - plt.legend() - plt.show() - # Analyze Protocol 1 π_minus_1 = 0.5 T_max = 30 @@ -1130,6 +1074,35 @@ $$ (eq:classerrorprob) where $\tilde \alpha_t = {\rm Prob}(l_t < 1 \mid f)$ and $\tilde \beta_t = {\rm Prob}(l_t \geq 1 \mid g)$. +Now let's write some code to simulate it + +```{code-cell} ipython3 +def compute_protocol_2_errors(π_minus_1, T_max, N_simulations, f_func, g_func, + F_params=(1, 1), G_params=(3, 1.2)): + """ + Compute error probabilities for Protocol 2. + """ + sequences, true_models = protocol_2(π_minus_1, + T_max, N_simulations, F_params, G_params) + l_ratios, _ = compute_likelihood_ratios(sequences, f_func, g_func) + + T_range = np.arange(1, T_max + 1) + + accuracy = np.empty(T_max) + for t in range(T_max): + predictions = (l_ratios[:, t] >= 1) + actual = true_models[:, t] + accuracy[t] = np.mean(predictions == actual) + + return { + 'T_range': T_range, + 'accuracy': accuracy, + 'l_ratios': l_ratios, + 'true_models': true_models + } + +``` + Since for each $t$, the decision boundary is the same, the decision boundary can be computed as ```{code-cell} ipython3 @@ -1177,11 +1150,11 @@ plt.tight_layout() plt.show() ``` -To the left of the green vertical line $g < f$, so $l_t < 1$; therefore a $w_t$ that falls to the left of the green line is classified as a type $g$ individual. +To the left of the green vertical line $g < f$, so $l_t > 1$; therefore a $w_t$ that falls to the left of the green line is classified as a type $f$ individual. - * The shaded orange area equals $\beta$ -- the probability of classifying someone as a type $g$ individual when it is really a type $f$ individual. + * The shaded red area equals $\beta$ -- the probability of classifying someone as a type $g$ individual when it is really a type $f$ individual. -To the right of the green vertical line $g > f$, so $l_t >1 $; therefore a $w_t$ that falls to the right of the green line is classified as a type $f$ individual. +To the right of the green vertical line $g > f$, so $l_t < 1$; therefore a $w_t$ that falls to the right of the green line is classified as a type $g$ individual. * The shaded blue area equals $\alpha$ -- the probability of classifying someone as a type $f$ when it is really a type $g$ individual. @@ -1213,6 +1186,29 @@ Now we simulate timing protocol 2 and compute the classification error probabili In the next cell, we also compare the theoretical classification accuracy to the empirical classification accuracy ```{code-cell} ipython3 +def analyze_protocol_2(π_minus_1, T_max, N_simulations, f_func, g_func, + theory_error=None, F_params=(1, 1), G_params=(3, 1.2)): + """Analyze Protocol 2.""" + result = compute_protocol_2_errors(π_minus_1, T_max, N_simulations, + f_func, g_func, F_params, G_params) + + # Plot results + plt.figure(figsize=(10, 6)) + plt.plot(result['T_range'], result['accuracy'], + 'b-', linewidth=2, label='empirical accuracy') + + if theory_error is not None: + plt.axhline(1 - theory_error, color='r', linestyle='--', + label=f'theoretical accuracy = {1 - theory_error:.4f}') + + plt.xlabel('$t$') + plt.ylabel('accuracy') + plt.legend() + plt.ylim(0.5, 1.0) + plt.show() + + return result + # Analyze Protocol 2 result_p2 = analyze_protocol_2(π_minus_1, T_max, N_simulations, f, g, theory_error, (F_a, F_b), (G_a, G_b)) @@ -1221,7 +1217,21 @@ result_p2 = analyze_protocol_2(π_minus_1, T_max, N_simulations, f, g, Let's watch decisions made by the two timing protocols as more and more observations accrue. ```{code-cell} ipython3 -# Compare both protocols +def compare_protocols(result1, result2): + """Compare results from both protocols.""" + plt.figure(figsize=(10, 6)) + + plt.plot(result1['T_range'], result1['error_prob'], linewidth=2, + label='Protocol 1 (Model Selection)') + plt.plot(result2['T_range'], 1 - result2['accuracy'], + linestyle='--', linewidth=2, + label='Protocol 2 (classification)') + + plt.xlabel('$T$') + plt.ylabel('error probability') + plt.legend() + plt.show() + compare_protocols(result_p1, result_p2) ``` @@ -1710,7 +1720,7 @@ P_g = np.array([[0.5, 0.3, 0.2], markov_results = analyze_markov_chains(P_f, P_g) ``` -## Related Lectures +## Related lectures Likelihood processes play an important role in Bayesian learning, as described in {doc}`likelihood_bayes` and as applied in {doc}`odu`. diff --git a/lectures/likelihood_ratio_process_2.md b/lectures/likelihood_ratio_process_2.md index fdff7ec2c..ba5c6a603 100644 --- a/lectures/likelihood_ratio_process_2.md +++ b/lectures/likelihood_ratio_process_2.md @@ -65,7 +65,7 @@ from math import gamma from scipy.integrate import quad ``` -## Review: Likelihood Ratio Processes +## Review: likelihood ratio processes We'll begin by reminding ourselves definitions and properties of likelihood ratio processes. @@ -166,7 +166,7 @@ def simulate(a, b, T=50, N=500): return l_arr ``` -## Blume and Easley's Setting +## Blume and Easley's setting Let the random variable $s_t \in (0,1)$ at time $t =0, 1, 2, \ldots$ be distributed according to the same Beta distribution with parameters $\theta = \{\theta_1, \theta_2\}$. @@ -195,7 +195,7 @@ $$c^1(s_t) = y_t^1 = s_t. $$ But in our model, agent 1 is not alone. -## Nature and Agents' Beliefs +## Nature and agents' beliefs Nature draws i.i.d. sequences $\{s_t\}_{t=0}^\infty$ from $\pi_t(s^t)$. @@ -239,7 +239,7 @@ $$ c_t^1 + c_t^2 = 1 . $$ -## A Socialist Risk-Sharing Arrangement +## A socialist risk-sharing arrangement In order to share risks, a benevolent social planner dictates a history-dependent consumption allocation that takes the form of a sequence of functions @@ -260,7 +260,7 @@ To design a socially optimal allocation, the social planner wants to know what a As for the endowment sequences, agent $i$ believes that nature draws i.i.d. sequences from joint densities $$ -\pi_t^i(s^t) = \pi(s_t)^i \pi^i(s_{t-1}) \cdots \pi^i(s_0) +\pi_t^i(s^t) = \pi^i(s_t) \pi^i(s_{t-1}) \cdots \pi^i(s_0) $$ As for attitudes toward bearing risks, agent $i$ has a one-period utility function @@ -269,7 +269,7 @@ $$ u(c_t^i) = \ln (c_t^i) $$ -with marginal utility of consumption in period $i$ +with marginal utility of consumption in period $t$ $$ u'(c_t^i) = \frac{1}{c_t^i} @@ -284,7 +284,7 @@ $$ (eq:objectiveagenti) where $\delta \in (0,1)$ is an intertemporal discount factor, and $u(\cdot)$ is a strictly increasing, concave one-period utility function. -## Social Planner's Allocation Problem +## Social planner's allocation problem The benevolent dictator has all the information it requires to choose a consumption allocation that maximizes the social welfare criterion @@ -303,8 +303,7 @@ This means that the social planner knows and respects * each agent's one period utility function $u(\cdot) = \ln(\cdot)$ * each agent $i$'s probability model $\{\pi_t^i(s^t)\}_{t=0}^\infty$ -Consequently, we anticipate that these objects will appear in the social planner's rule for allocating the aggregate endowment each period. - +Consequently, we anticipate that these objects will appear in the social planner's rule for allocating the aggregate endowment each period. First-order necessary conditions for maximizing welfare criterion {eq}`eq:welfareW` subject to the feasibility constraint {eq}`eq:feasibility` are @@ -312,9 +311,6 @@ $$\frac{\pi_t^2(s^t)}{\pi_t^1(s^t)} \frac{(1/c_t^2(s^t))}{(1/c_t^1(s^t))} = \fra which can be rearranged to become - - - $$ \frac{c_t^1(s^t)}{c_t^2(s^t)} = \frac{\lambda}{1-\lambda} l_t(s^t) $$ (eq:allocationrule0) @@ -361,7 +357,7 @@ $$ -## If You're So Smart, $\ldots$ +## If you're so smart, $\ldots$ Let's compute some values of limiting allocations {eq}`eq:allocationrule1` for some interesting possible limiting @@ -400,7 +396,7 @@ Doing this will allow us to connect our analysis with an argument of {cite}`alch -## Competitive Equilibrium Prices +## Competitive equilibrium prices Two fundamental welfare theorems for general equilibrium models lead us to anticipate that there is a connection between the allocation that solves the social planning problem we have been studying and the allocation in a **competitive equilibrium** with complete markets in history-contingent commodities. @@ -526,6 +522,9 @@ According to formula {eq}`eq:pformulafinal`, we have the following possible limi * when $l_\infty = \infty$, $c_\infty^1 = 1 $ and tails of competitive equilibrium prices reflect agent $1$'s probability model $\pi_t^1(s^t)$ according to $p_t(s^t) \propto \delta^t \pi_t^1(s^t) $ * for small $t$'s, competitive equilibrium prices reflect both agents' probability models. +We leave the verification of the shadow prices to the reader since it follows from +the same reasoning. + ## Simulations Now let's implement some simulations when agent $1$ believes marginal density @@ -821,7 +820,7 @@ This ties in nicely with {eq}`eq:kl_likelihood_link`. -## Related Lectures +## Related lectures Complete markets models with homogeneous beliefs, a kind often used in macroeconomics and finance, are studied in this quantecon lecture {doc}`ge_arrow`. @@ -833,11 +832,6 @@ Likelihood ratio processes appear again in {doc}`advanced:additive_functionals`. -{doc}`ge_arrow` - - - - ## Exercises ```{exercise} diff --git a/lectures/likelihood_var.md b/lectures/likelihood_var.md index 2bacf4812..9f76cf2b0 100644 --- a/lectures/likelihood_var.md +++ b/lectures/likelihood_var.md @@ -20,7 +20,7 @@ kernelspec: ``` -# Likelihood Processes for VAR Models +# Likelihood Processes For VAR Models ```{contents} Contents :depth: 2 @@ -156,11 +156,11 @@ Given the Gaussian structure, the conditional distribution $f(x_{t+1} | x_t)$ is - Mean: $A x_t$ - Covariance: $CC'$ -The log conditional density is: +The log conditional density is $$ \log f(x_{t+1} | x_t) = -\frac{n}{2} \log(2\pi) - \frac{1}{2} \log \det(CC') - \frac{1}{2} (x_{t+1} - A x_t)' (CC')^{-1} (x_{t+1} - A x_t) -$$ +$$ (eq:cond_den) ```{code-cell} ipython3 def log_likelihood_transition(x_next, x_curr, model): @@ -234,7 +234,26 @@ def simulate_var(model, T, N_paths=1): ## Likelihood ratio process -Now let's compute likelihood ratio processes for comparing two VAR models +Now let's compute likelihood ratio processes for comparing two VAR models. + +For a VAR model with state vector $x_t$, the log likelihood ratio at time $t$ is + +$$ +\ell_t = \log \frac{p_f(x_t | x_{t-1})}{p_g(x_t | x_{t-1})} +$$ + +where $p_f$ and $p_g$ are the conditional densities under models $f$ and $g$ respectively. + +The cumulative log likelihood ratio process is + +$$ +L_t = \sum_{s=1}^{t} \ell_s = \sum_{s=1}^{t} \log \frac{p_f(x_s | x_{s-1})}{p_g(x_s | x_{s-1})} +$$ + +where $p_f(x_t | x_{t-1})$ and $p_g(x_t | x_{t-1})$ are given by their respective conditional densities defined in {eq}`eq:cond_den`. + + +Let's write those equations in Python ```{code-cell} ipython3 def compute_likelihood_ratio_var(paths, model_f, model_g): @@ -271,7 +290,7 @@ def compute_likelihood_ratio_var(paths, model_f, model_g): return log_L_ratios if N_paths > 1 else log_L_ratios[0] ``` -## Example 1: Two AR(1) processes +## Example 1: two AR(1) processes Let's start with a simple example comparing two univariate AR(1) processes with $A_f = 0.8$, $A_g = 0.5$, and $C_f = 0.3$, $C_g = 0.4$ @@ -314,7 +333,7 @@ plt.show() As we expected, the likelihood ratio processes goes to $+\infty$ as $T$ increases, indicating that model $f$ is chosen correctly by our algorithm. -## Example 2: Bivariate VAR models +## Example 2: bivariate VAR models Now let's consider an example with bivariate VAR models with