May 26, 2026

[DRAFT] Two Probabilistic Models of NIR-2025

This post introduces probabilistic modelling to PlanZero. Probabilistic modelling extends simulation-based modelling with a capacity to fit model parameters to data. Two models are introduced: a simple baseline that assumes emissions don't change over time, and a more sophisticated model based on a two-step autoregressive process (AR-2). Both models estimate provincial and territorial emissions for every sector and greenhouse gas, in which provincial and territorial emissions sum up to national totals.

Table of Contents:

Introduction
Prediction with Quantified Confidence: the Motivation for Probabilistic Modelling
Probabilistic Modelling of Harvested Wood Products' ${C O}_{2}$ with Static Distributions
Modelling of Harvested Wood Products ${C O}_{2}$ with an Autoregressive Process
- Inference of Harvested Wood Products ${C O}_{2}$ emissions with an AR2-based model
- Inference of All Emissions with an AR2-based model
Next Steps and Conclusion

Introduction

If I were to predict that the Canadian national emissions total for 2025 would be exactly 690.395 Mt ${C O}_{2} e$ , it would be silly because (a) how could the total be that exact, specific number, and (b) who would care about that degree of precision when there are a few percentage points worth of error on those totals anyway?

If I were to claim instead that they might be exactly that amount, then again it would be a silly statement because of course it's possible; why not?

The way people have learned to make predictions about such things is to use the language of probability. A probability distribution is a table (or function) of the relative probabilities of all the possible values of an unknown quantity (a random variable), such as the 2025 emissions total in the future 2027 National Inventory Report. The emissions total has to be some value so the probabilities of all of it's possible values add up to 1. A non-silly prediction about the 2025 emissions in the 2027 report might be that "there's a 15% chance the total will be between 690 and 700 Mt ${C O}_{2} e$ ." That prediction would be the result of adding up the probabilities of values from 690 to 700 according to a probability distribution.

PlanZero has not, until this post, produced probability distributions for e.g. future emission quantities. The post introduces probabilistic modelling via two models: a trivial baseline model that ignores any temporal trends in the data whatsoever (the static model), and non-trivial baseline model based on a two-step autoregressive process that can match trends (the AR-2 model). The static model is simpler, and serves to explain how probabilistic modelling works. Its random variables are easier to visualize. The AR2 model is more complex, but provides a better fit to the data.

Probabilistic Modelling of National Emissions with Static Normals

My mental model of the 2025 National Inventory Report (NIR-2025, or just NIR in this post), is that for each

year, from 1990 to 2023,
province and territory (I'll call these "regions"),
IPCC sector (I'll call these just "sectors"),
and greenhouse gas,

there is a quantity, in units of e.g. kt

{C O}_{2} e

that has been estimated in support of Canada's reporting obligations under the Paris Agreement. Let's call this "hypothesis 1". Furtheremore, I believe ("hypothesis 2") that the national totals in the NIR are sums of subtotals from provinces and territories (regions).

Alas, neither of these hypotheses is strictly true. They are true for the majority of years, territories, sectors, and gases, but there are counter-examples to both. Counter to hypothesis 1, there are missing records in the data, such as 81 relating to stationary combustion in the construction industry. Counter to hypothesis 2, regional totals do not always sum up to the national one. If we define the sum of regions vs. national total discrepancy as the maximum difference over years and provinces, there are 11 with a discrepancy over 1Mt ${C O}_{2} e$ , which are listed in the table below (out of 486 sector and GHG combinations). At the level of 100kt ${C O}_{2} e$ discrepancy, there are about twice as many records, and at the level of 10kt ${C O}_{2} e$ discrepancy, twice as many again. I defined "Relative Discrepancy" as the discrepancy divided by the greater of itself and the root mean square of national totals, as a way to relate discrepancies to the scale of numbers in their sectors.

IPCC Sector	GHG	Max Discrepancy (kt ${C O}_{2} e$ )	Relative Discrepancy
Petroleum Refining Industries	${C O}_{2}$	2430.85	14.4%
Commercial and Institutional	${C O}_{2}$	1271.44	4.1%
Cement Production	${C O}_{2}$	2421.59	36.7%
Lime Production	${C O}_{2}$	1421.05	84.4%
Ammonia Production	${C O}_{2}$	3090.34	100.0%
Nitric Acid Production	$N_{2} O$	1097.89	100.0%
Petrochemical and Carbon Black Production	${C O}_{2}$	4398.65	100.0%
Non-Energy Products from Fuels and Solvent Use	${C O}_{2}$	7395.02	67.7%
Non-Energy Products from Fuels and Solvent Use	$N_{2} O$	1068.83	100.0%
Forest Land	${C O}_{2}$	1731.42	2.1%
Harvested Wood Products	${C O}_{2}$	53296.13	100.0%

For several modest-to-minor sectors, the size of the discrepancy is sometimes at the same order of magnitude as the national total, and even for large sectors such as harvested wood products and [static combustion in] commercial and institutional environments, there are significant discrepancies.

Inference in probabilistic models is about finding values for unobserved random variables that balance two competing objectives: to find values that are probable according to the prior beliefs which are also consisent with observation. The data in the NIR-2025 is not meant to be interpreted too literally. In fact, NIR-2025 includes per-sector estimates of uncertainty on national totals for the base year (1990) and the most-recent year (2023). These uncertainties should be considered a lower bound on our uncertainty about the true sectoral emissions, because for example it can also happen, from time to time, that the protocols for estimating emissions are updated. This happened a few years ago in the estimation of vented ${C H}_{4}$ , which caused a re-assessment of the past estimates, and a change that went beyond the previously-estimated uncertainty.

So if we're missing data, and the data we have in the NIR doesn't add up, and anyway the data is published with lower-bounded uncertainty, then what can we say about what really happened? To answer this question, we have to predict what probably happened, and we do this by building a probabilistic model. Such a model consists of

random variables for what we observe (e.g. data in the NIR),
random variables for what we think should be there (e.g. neat, coherent emissions totals that all add up etc.),
mathematical relationships between them (e.g. NIR values might be Normally-distributed around the neat, coherent values).

Within the context of a given model, we can infer the likely values of the neat, coherent variables given what we observe in the NIR data, automatically, using rules of probability theory. (Whether the model is actually good, or which of two models is better, is an important question that I'll come back to in a future post, but leave un-addressed for now.) The next section defines what I'm calling the static model, beginning with a brief refresher of probability notation.

Modelling Harvested Wood Products with Static Normals

One greatly-simplifying modelling hypothesis we could make is that sectoral emissions (per territory, per GHG) aren't actually changing at all over time. With this assumption we can go through the steps of modelling and analysis for a model that is relatively easy to illustrate, compute, etc. Later we will compare it with a model that recognizes temporal trends. This section expresses a probabilistic model of emissions.

Some background on probability notation is required, in order to express probability models.

$P (X = 0)$ means there's a number we don't know, we're calling it $X$ , and this is the probability it is actually 0. $X$ is called a random variable.
$P (X)$ is the table (possibly infinitely large, in which case it might be called a function) of all the values $X$ might take, and their probabilities (or densities, if the set of values $X$ might take is continuous)
$P (0)$ is sloppy shorthand for $P (X = 0)$ when it's clear from surrounding context we're talking about $X$
$𝒩 (μ, σ^{2})$ is a Normal distribution, which is non-negative function that sums (integrates) to 1. It is parameterized by two numbers: $μ$ the mean, and $σ^{2}$ the variance. (Sometimes these are called Gaussians but I'm going to call them Normals to be consistent with the probabilistic programming libraries I'm aware of.)
$X ~ 𝒩 (0, 1)$ means $X$ is a random variable, and its corresponding distribution $P (X)$ is the thing to the right of the " $~$ " symbol, which is a mapping from values to probabilities (or densities if the values are continuous instead of discrete).
$P (A = 0, B = 1)$ means there's multiple numbers we don't know, and this is the probability that they both actually have these specific values.
$P (A, B) = P (A) P (B)$ means that the probability distributions for $A$ and $B$ are independent. A distribution over multiple random variables is called a joint distribution whether they are independent or not.
$P (A = 0 | B = 1)$ is the probability that $A = 0$ if we make an additional assumption about $B$ , namely that it is 1. If A and B are not independent, then knowing something about B tells us something about A; in other words, changes what would otherwise have been our $P (A)$ .

Probabilistic modelling is a huge topic. I think it's typically introduced in secondary school, and covered in greater detail at various university levels, often with application to areas of science, engineering, medicine, social studies, mathematics etc. where probabilities are used. For more information and orientation, see e.g. wikipedia ( Statistical Model, Graphical Models), many tutorials online specific to probabilistic programming software packages (e.g. Stan, PyMC, Pyro; PlanZero is currently using NumPyro but that might change in the future) or a textbook such as Probabilistic Graphical Models (Koller and Friedman, 2009) for greater depth.

With the two-paragraph crash course in probability out of the way, the following table expresses the static (non-time-varying) probabilistic model model that will be our subject for the rest of this section:

$μ_{(s, g, r)} ~ 𝒩 (0, 1)$	Suppose the regional mean emissions are probably around 0, plus or minus 1. It is customary to scale data to be in this range before modelling them, so each emissions time series is scaled to have a mean squared value of 1. After this scaling, plus-or-minus 1 is a relatively un-opinionated prior as to what $μ$ should be, relative to the NIR data, which will probably say that $μ$ should be somewhere in the range of plus-or-minus 1.
$σ_{(s, g, r)} ~ L o g N o r m a l (- 1, 0.7)$	Suppose regional variances about the mean are never super-tiny, nor much larger than 1.
$σ_{(s, g)} ~ L o g N o r m a l (- 1, 0.7)$	Suppose national variance about the mean are never super-tiny, nor much larger than 1.
$X_{(s, g, r, t)} ~ 𝒩 (μ_{(s, g, r)}, σ_{(s, g, r)}^{2})$	Suppose that for each sector $s$ , greenhouse gas $g$ , region $r$ (province or territory), and time $t$ (annual total), suppose there's a number we don't know (random variable X) that's normally distributed about some regional mean $μ_{(s, g, r)}$ , and with some variance $σ_{(s, g, r, t)}^{2}$ .
$Y_{(s, g, t)} ~ 𝒩 (\sum_{r \in r e g i o n s} μ_{(s, g, r)}, σ_{(s, g)}^{2})$	Suppose further that for each sector $s$ , greenhouse gas $g$ , and time $t$ (annual total), some number we don't know (random variable Y) is normally distributed about a mean that is the sum of the regional means, with a variance $σ_{(s, g)}^{2}$ .

The prior densities for the

μ

and

σ

terms are illustrated below.

That's it, that's a model of an NIR! The next section is about how we can use an inference algorithm and NIR-2025 data to predict likely neat and coherent latent values ( $μ$ ), albeit in the context of this static model.

Inference Algorithms and Probabilistic Programming

The prior model is not really meant to be very accurate, even considering its structural limitation of being a static model. To fix it up, the expectation is to perform some amount of inference. The natural inference for this model is to treat the actual NIR as observations of our $X$ and $Y$ random variables and infer what $μ$ and $σ$ probably were, given both our prior about what they would be, and the observation of X and Y (which depend on $μ$ and $σ$ ). Bayesian inference addresses the question of what $μ$ and $σ$ probably were in terms of a joint probabilility distribution:

$P (μ, σ | X = N I R^{(2025, r e g i o n a l)}, Y = N I R^{(2025, n a t i o n a l)})$

This [joint] distribution is specifically a posterior joint distribution because it's conditioned on data. Latent [random] variables are the variables in our model that are not observed, i.e. for which we don't have data. Generally, the posterior distribution over a model's latent variables is not a distribution with a neat and tidy formula such as a normal or exponential. Typically the only practical way to work with these distributions is with algorithms, such as the inference algorithms mentioned in the previous section. Inference algorithms work on the basis of Bayes' Theorem. In standard probability notation, Bayes' theorem looks like:

$P (A | B) = \frac{P (A, B)}{P (B)} = \frac{P (B | A) P (A)}{P (B)}$

Bayes' theorem shows how to relate distributions of the form $P (A | B)$ , to ones involving $P (B | A)$ instead. With the help of Bayes theorem, it is possible to transform the posterior distribution we're interested in, $P (μ, σ | X = N I R^{(2025, r e g i o n a l)}, Y = N I R^{(2025, n a t i o n a l)})$ , into a formula with the sorts of pieces with which we specified our model:

$\frac{P (X = N I R_{(2025, r e g i o n a l)}, Y = N I R_{(2025, n a t i o n a l)} | μ, σ) P (μ) P (σ)}{P (X = N I R_{(2025, r e g i o n a l)}, Y = N I R_{(2025, n a t i o n a l)})}$

Probabilistic programming systems enable someone such as myself to write down a generative sampling processes (such as the one in the left column of the table above that defined the model) in a programming language, and leverage standard algorithm implementations to infer samples from the posterior distribution. Even with a Ph.D. in machine learning, I still find these algorithms amazing (read up on e.g. MCMC algorithms or stochastic variational inference if you're interested in how these algorithms work). The implementation of these algorithms can involve relatively sophisticated software techniques (e.g. tracing and runtime program transformations) and the algorithms themselves benefit from ongoing statistics research. For our purposes we'll take their capabilities as given (specifically via the use of Jax and NumPyro, and the NUTS MCMC algorithm), and apply them, in the next section, to draw samples from the posterior distribution $P (μ, σ | X = N I R^{(2025, r e g i o n a l)}, Y = N I R^{(2025, n a t i o n a l)})$ assuming the static model described above.

Inference of Latent Variables in the Static Model

Running NumPyro's NUTS inference algorithm for the static model above, using the NIR-2025 emissions data regarding ${C O}_{2}$ emissions from the Harvested Wood Products sector (HWP), results in the following ranges of values for $μ_{(H W P, C O_{2})}$ (shown as darker shaded region) and $σ_{(H W P, C O_{2})}$ (shown via lightly-shaded 95% credible intervals on the regional and per-year random variables within $X_{(H W P, C O_{2})}$ and $Y_{(H W P, C O_{2})}$ ).

Looking at the modelling results ( $μ$ and $σ$ ) we see some good and some bad. It is good, that $μ$ (mean is dotted horizontal line, 95% credible interval is darker shaded region) typically goes through the data points somewhere near the middle. Likewise, $σ$ has been estimated such that the data mostly fall within the lightly-shaded regions, with a few points popping out above or below; which is approximately we expect from a 95% credible interval after fitting to the data. That said, it is bad that we see large blocks of shaded region without any data, which is a symptom of temporal correlation in the data. (We'll compare this quality of fit with a different model in the next sections.)

Note also that there are differences between the horizontal solid lines (each region's sample mean) and the horizontal dotted lines (each region $r$ 's mean $μ_{(H W P, C O_{2}, r)}$ ). The darker shaded region helps put their difference into perspective. The darker shaded region is the 95% credible interval for $μ_{(H W P, C O_{2}, r))}$ (the upper and lower limits of the central 95% of sampled values from the latent variable's posterior distribution). For two regions, British Columbia and Quebec, the data mean lies slightly outside of their regional $μ_{(H W P, C O_{2}, r)}$ credibility interval. This means that it is improbable, based on the model, and the data, that either $μ_{(H W P, C O_{2}, B C)}$ or $μ_{(H W P, C O_{2}, Q C)}$ is really the data mean for each of those provinces. I chose Harvested Wood Products' ${C O}_{2}$ emissions because there were a number of years where the regional and national totals are thousands of kt ${C O}_{2}$ apart. Unsurprisingly, inference in this model reveals that a discrepancy also exists terms of sample means-since-1990. Remember, the regional latent random variables $μ_{s, g, r}$ in the model were required to add up to a national total, whereas actual regional and national totals in NIR-2025 do not always add up this way. The inference algorithm has estimated a sample of values for $μ_{s, g, r}$ that balances (a) fitting the regional emissions data, (b) fitting the national total emissions, and (c) respecting the prior distributions in the model.

That concludes PlanZero's first probabilistic model of a sector and gas (Harvested Wood Products ${C O}_{2}$ ). A standard view of probabilistic models of sectors now features this Static Normals model for Harvested Wood Products on a page that shows the following figure:

This figure doesn't show $μ$ specifically because it's a standard style of figure for any of probabilistic model of that sector. But it does show the posterior samples of $μ$ $σ$ relate to the set of plausible emissions values.

In fact, PlanZero now includes this sort of modelling for all of Canada's emissions sectors. The Static Normals model overview includes the following figure, showing some quality of fit across the board. The colours correspond to high-level emissions categories:

dark blue for energy used in stationary combustion (e.g. for public electricity, oil and gas extraction, etc.)
light green for energy used in vehicles (e.g. light-duty, heavy-duty, off-road, marine, etc.)
purple for fugitive sources(e.g. venting, flaring)
wheat for agriculture-related emissions (e.g. soil, manure)
pink for waste-related emissions (e.g. landfills)
other colours for various industrial activities (mineral products, chemicals, metals)

The emissions sectors are ordered according to decreasing upper bound on their credible intervals. The largest sectors (e.g. public heat, oil and gas extraction) are across the top after the national total without Land-Use, Land-Use Changes, and Forestry (LULUCF), on a common scale from 0 to 150 Mt

{C O}_{2} e

. The next row down has emissions on a scale from 0 to 52 Mt

{C O}_{2} e

with things like residential heating and enteric fermentation, and subsequent rows go down from there to the bottom row where the scale is relatively small.

Mouse-over the slider below to see the active date range and adjust it to focus on specific periods. Click on each panel for a deeper-dive into modelled emissions in that sector.

That concludes the introduction of PlanZero's "Static Normals" model. It is a probabilistic model, but the modelling assumption that each NIR emission sector total is actually unchanging throughout history and into the future is a very untrue assumption for many sectors. There has been consistent significant reduction in emissions from e.g. public electricity production over the past two decades as coal generators have been phased out. On the other hand, the emissions from oil and gas extraction have risen enormously over that same period due to changes in the scale and methods of oil extraction operations. The Static Normals model is designed for simplicity, and it cannot model such effects. The next section introduces a different model based on auto-regressive random processes, which has the property that inference can be described as learning trends in data (this is why machine learning in the field of computer science has that name). The predictions of the AR2 model are more aligned with at least my own, if not most people's intuitions, for how emissions rates change over time.

An Alternative Model based on an Autoregressive Process

It is natural to view the emissions data in the NIR as time series with trends. A model should leverage the fact that emissions tend to change slowly; the emissions from the year before are usually a good guess for the current year, and probably a better guess than the overall mean since 1990 (which is what the Static Normals model does). There are many approaches and applications of time-series analysis but this section develops an emissions model based on autoregressive analysis. An auto-regressive process of order $p$ is a set of random variables that can be written like this:

$\begin{matrix} Z_{1 . . p} & ~ & 𝒩 (0, σ^{2}) \\ Z_{t} & ~ & 𝒩 (c_{0} + c_{1} Z_{t - 1} + c_{2} Z_{t - 2} + . . . + c_{p} Z_{t - p}, σ^{2}) \end{matrix}$

There's an order to them. The first few random variables are just independent normals (the first $p$ of them: $Z_{1 . . p} ~ 𝒩 (0, σ^{2})$ ) but after that each one is a weighted combination of the previous ones plus a bit of noise. Sometimes this sort of set of random variables is called a random walk, because each one takes a random step along a trend that's established by the weights $c$ . Whether the random walk is dominated by the trend or the noise depends on the magnitude of the weights $c$ and the magnitude of the noise $σ$ .

The following time-varying graphical model for NIR data will be our subject for the rest of this section. It is based on $μ_{(s, g, r, t)}$ being a two-step auto-regressive process. Two-step auto-regressive processes can model damped systems with momentum (I chose two steps because one step seemed too simple — just damping.) Indexes by $s, g$ are dropped in the table below for brevity, and because anyway all sectors and gases are considered to be independent. So for example what is actually $μ_{(s, g, r, t)}$ in the context of the national model is written just $μ_{(r, t)}$ .

$c_{(r, 1)} = 3 K u m a r a s w a m y (1.25, 1.25) - 0.5$	Let the coefficient on the previous timestep range between -0.5 and 2.5. The Kumaraswamy distribution with parameters (1.25, 1.25) provides a distribution that's a bit more uniform-looking than the (2,2) distribution illustrated on the wikipedia page. Parameters (2,2) may have worked here too, I didn't try them.
$c_{(r, 2)} = 3 K u m a r a s w a m y (1.25, 1.25) - 1.5$	Let the coefficient on the 2nd-previous timestep range between -1.5 and 1.5.
$σ_{r} = max (\sqrt{\frac{1}{T} \sum_{t}^{T} e_{r, t}^{2}}, 1)$	Even after scaling the sector emissions data set so that the national totals had a dynamic range of about 1, still peek at the provincial and territorial emissions data to define a scale for each of these regions, because some provinces are way more active in some sectors than others. (The "e" is for the regional emissions after scaling by the national factor.)
$\forall_{r} \forall_{t \in {1, 2}} μ_{(r, t)} ~ 𝒩 (0, σ_{r}^{2})$	Let the first two steps of the autoregressive process for each region be independent normals, scaled appropriately for the region.
$μ_{(r, t)} ~ 𝒩 (c_{(r, 1)} μ_{(r, t - 1)} + c_{(r, 2)} μ_{(r, t - 2)}, {0.05}^{2} σ_{r}^{2})$	Suppose there's a sequence of random variables indexed by time $μ_{(r, t)}$ that unfold through time according to the AR-2 process introduced above. The ${.05}^{2} σ_{r}^{2}$ term here controls how quickly the AR-2 process can change course from where its previous values and weights $c$ would otherwise take it. Inference is quite sensitive to this term.
$X_{(r, t)} ~ S B L N (μ_{(r, t)}, {0.05}^{2} σ_{(r, t)}^{2})$	Suppose there are random variables indexed by region and year $t$ that have their own variance, but whose means are identical to $μ_{(r, t)}$ . These are the variables that are meant to represent the per-region (and per-gas, per-sector) emissions totals in the NIR, whenever the NIR includes such totals. The variance in these model was fixed, and based roughly on the uncertainty estimate in the NIR-2025. The SBLN distribution, described below, is like a Normal distribution when $μ$ is around 0, more like a Log-Normal distribution when $μ$ is large and positive, and more like a negative Log-Normal when $μ$ is large and negative.
$Y_{t} ~ S B L N (\sum_{r \in r e g i o n s} μ_{(r, t)}, {0.05}^{2})$	Finally, suppose that there are time-indexed random variables with their own variance, whose means are the sum of regional variables $μ_{(r, t)}$ , which represent the (per-gas, per-sector) national totals. Similarly to $X$ , variance here was fixed, based roughly on the uncertainty estimate for NIR-2025.

The $S B L N (μ, σ^{2})$ distribution is introduced here as a hybrid between a Normal and a Log-Normal distribution. Uncertainty in the NIR is typically proportional to the size of the emissions amounts in question (which vary widely in magnitude), and is reported in terms of percentage points, which suggests a Log-Normal distribution for e.g. $X | μ$ might be more appropriate than a Normal. Except that sometimes, especially in Harvested Wood Products, the emissions are sometimes negative and sometimes positive, and a Log-Normal distribution can only do negative or positive. To get around this, I introduced a special distribution that I call a Symmetric Blended Log-Normal that gradually transitions from being a negative Log-Normal, to a regular Normal, to a positive Log-Normal as its $μ$ parameter sweeps from being large and negative, through zero, and to being large and positive. An implementation and illustration of the distribution in a Jupyter notebook is rendered nicely on GitHub, for anyone looking for more detail about this.

Inference of Harvested Wood Products ${C O}_{2}$ emissions with an AR2-based model

Running NumPyro's NUTS inference algorithm for the AR2 model above, using the NIR-2025 emissions data regarding ${C O}_{2}$ emissions from the Harvested Wood Products sector (HWP), results in the following values for $μ_{(H W P, C O_{2})}$ (shown via dotted mean and shaded region credible intervals, scaled to match the data magnitude).

In contrast to the Static Normals model, the AR2 model predicts different emissions at different times. Around the data points, the credibility intervals usually lie near the data, and are pretty narrow, but don't always include the data. HWP is not an easy running example, but it is an interesting one. The assumptions of the AR2 model are in conflict with the reality of HWP emissions. HWP emissions jump around a lot rather than being smooth, and the national total is almost 30 Mt ${C O}_{2} e$ off the sum of regional subtotals in the 1990s. Still, the goal of inference is to find the most-plausible values for the latent random variables, especially $μ$ , and as far as I can tell the NUTS algorithm is working.

The AR2 model is a better starting point for future work than the Static Normals model, but more work is in order. For example, here in the real world beyond the model, it is unlikely that the Northwest Territories and Nunavut were secretly emitting vast quantities of HWP ${C O}_{2}$ in the early 90s, in the years where no emissions data are listed for them in NIR-2025. Inference has produced an interesting theory that helps explain the regional sum vs. national total discrepancy, but it is not true. (Nunavut has a good excuse for not listing emissions: it didn't exist until 1999.) For this reason and others, the AR2-based model of HWP emissions could be improved by e.g.

stronger priors regarding historical emissions when NIR data is missing
incorporation of the mechanisms of HWP calculations (which work in terms of filters on quantities of types of harvested wood, see Emissions Calculations for Harvested Wood Products and Forest Land (Feb 22, 2026) for an earlier review)
combined modelling of both Forest Land and HWP emissions
the inclusion of more data sources regarding what was harvested, when, and for what purposes
other changes reflective of why the provincial and national estimates are so different

I expect to get to some of these points as I work back over IPCC sectors studied earlier this spring, and incorporate data sources and learnings from each one into future iterations of the AR2 probabilistic model.

For now, PlanZero's standard view of probabilistic models by sector includes this AR2 model for Harvested Wood Products figure:

It's fun to play with the slider to see what's going on in the 90s, and the panel view highlights that future uncertainty in the AR2 model is dominated by BC, ON, QC, and AB.

Inference of All Emissions with an AR2-based model

The type of figure used above to visualize predictions of Canada's future emissions according to the Static Normals model can also be used for the AR2 model. Again the sectors are ordered according to decreasing upper bound on their credible intervals, but because the models make different predictions, the ordering of sectors is slightly different. The figure is less colourful as well, because the colours highlight uncertainty, and the AR2 model predicts tighter, smaller, 95% credible intervals. Feel free to click any of the panels for a closer look at regional and per-gas emissions inferences and predictions. This figure is featured on the PlanZero AR2 overview page.

Overall, I think the AR2 predictions are reasonable. In most cases, they look like visual extensions of apparent trends, which is the sort of thing an AR2 model can be expected to do. There are some examples of poor fit, such as Ammonia Production, Nitric Acid Production, and Petrochemical and Carbon Black Production where the inferred total is far from the national data, but these cases are a bit like the case of Harvested Wood Products, in which the provincial and territorial sums didn't match the national total. In fact, in some of these cases, the provincial emissions are listed as all being zero. Inference in the AR2 model understandably struggles to explain how all of the regional emissions related to Ammonia Production can be zero while the national total is about 2.5 Mt ${C O}_{2} e$ . I hope modelling in these sectors can be improved by a better understanding on my part of what's going on, I haven't studied these sectors at all yet.

Next Steps and Conclusion

This introduction of probabilistic modelling to PlanZero was motivated by three things:

Potential PlanZero users have expressed interest in short-term forecasting.
Writing the glossary post about the computational and semantic aspects of modelling without mention of probabilistic modelling felt wrong.
I felt that predictions of the kind attempted in the Bovaer post would not be compelling without acknowledging and quantifying uncertainty, especially for strategies that are more active.

This post makes meaningful progress on all three fronts. The AR2 model is a full predictive model of the entire NIR. PlanZero's mission called for such a thing, and it didn't have one before. This full predictive model can be used for short-term, medium-term, and long-term forecasting. The AR2 model presented here leaves room for improvement in forecasting accuracy, but it provides a non-trivial baseline upon which to improve. The glossary post (still in Draft) can now be completed with a full overview of probabilistic and non-probabilistic modelling. Finally, future strategy analysis can build on the AR2 probabilistic model. This is especially valuable for active strategies, for which there is a little data and therefore an opportunity for posterior inference. Questions such as "How efficient is it to subsidize EV purchases?" for example, is a question that is better-suited to probabilistic modelling with Bayesian inference than to the non-probabilistic kind used to analyze the efficiency of Bovaer, because our prior beliefs about decision-making across the market can be refined by month-by-month data on uptake.

The PlanZero site has grown, with this post, to include a "Models" tab, which replaces the old "Simulations" tab, and moves the simulations themselves into a new section on that page for "Simulation-based models". That page also has a section for "Probabilistic Models" with brief descriptions and links to the Static Normals and AR2 models introduced in this post. The models predict emissions out to 2050 on the basis of 95% credibility intervals, with the Static Normals model predicting 379.2 - 1106.7 Mt ${C O}_{2} e$ in 2050 and the AR2 model predicting 737.7 - 862.9 Mt ${C O}_{2} e$ .

This post has also brought several new terms to the Glossary from foundational terms such as Probability Distribution and Random Variable to more advanced ones such as Probabilistic Model, Bayesian Inference, and Credible Interval.

Next Steps

The introduction of probabilistic modelling to PlanZero opens up several possible directions for next steps: inference diagnostics, prediction challenges and model comparison, model improvements.

Inference Diagnostics

The uncertainty estimates on latent variables in probabilistic models are not guaranteed to be good estimates. The algorithms that make these estimates can fail, and there are techniques for monitoring and communicating how well they've worked. I've made a GitHub issue to show such metrics alongside standard probabilistic modelling results after a few upgrades to AR2. It may be natural to extend the notion of Barrier to probabilistic models, and use per-barrier visualizations as the framing for inference diagnostics. Correct inference is more important when predictions have to be good, so this next step can wait until after PlanZero has established some prediction challenges, as described in the next section.

Prediction Challenges and Model Comparison

Probabilistic models can be used for a multitude of purposes. Typically many models are candidates for any given purpose. The question of which model to use for which purpose is not always easy to answer. The field of machine learning has made good progress developing better and better models for decades, by aligning around a set of standard prediction challenges (each with specific model-scoring systems) and iteratively improving performance on the standard challenges. One direction for future work on PlanZero is to add standard prediction challenges. In fact, I think this is the most important next step, because it is a necessary foundation for the evaluation of the various model improvements suggested in the next section.

Model Improvements

Even without model-evaluation metrics at hand, I'm going to go ahead and say that the AR2 model presented here leaves room for improvement (even if it's already better than the Static Normals). For each of the sectors I looked at this spring, its AR2 model can be extended to predict not just the NIR emissions, but also the data I found to explain the NIR emissions; the more types of data that $μ$ can explain via a probabilistic model, the tighter an estimate we can make for what $μ$ really was at various times up to the present, and how it should evolve forward in time. There are structural dynamics in the relationship between Harvested Wood Products emissions and Forest Land data that haven't yet been used to inform the AR2 model's predictions. There are natural dynamics in the life cycles of cattle that can relate StatsCan's cattle survey data to future methane emissions. There are natural dynamics in vehicle ownership cycles that can relate StatsCan's vehicle purchase data to future ${C O}_{2}$ emissions. Climate change and population change forecasts can help explain trends in several sectors in a way that can be extrapolated with greater confidence into the future.

Probabilistic modelling provides a framework within which all of these factors and mechanisms can be incorporated and analyzed together or individually as appropriate. There are many possible next steps that may be motivated as improvements to PlanZero's probabilistic models. My next step will be to implement prediction challenges that can serve as scaffolding for the iterative development of useful models.

Conclusion

Probabilistic modelling sets a new standard for what it means to add NIR sectors to PlanZero, as posts have done throughout much of Spring 2026. Probabilistic models together with Bayesian inference can be more credible and engaging than simulations because they can start with people's existing beliefs, data they can recognize as valid, and then proceed automatically to conclusions that they may not have considered. The probabilistic models added to PlanZero leave much room for improvement, but are already interesting to explore. The next post will introduce a prediction challenge that will serve to rank the two models here (Static Normals vs. AR2) as predictors of upcoming NIR results, and set a bar for future probabilistic model improvement.

Until then,

- James Bergstra