Conditional Likelihood

The ConditionalLikelihood algorithm computes the prediction error decomposition log-likelihood for fully-observed state-space models. At each time step, it predicts the next observation from the observed current state (not the model-predicted state), and accumulates the Gaussian log-likelihood of the innovation (prediction error).

This is the standard approach for maximum likelihood estimation of AR(1), VAR(1), nonlinear DSGE, and other models where the state is directly observed.

When to Use Each Algorithm

Mathematical Formulation

Given a state-space model with transition $x_{t+1} = f(x_t, w_t)$ and observation $z_t = g(x_t)$, the conditional log-likelihood is:

\[\log L = \sum_{t=1}^{T} \left[ -\frac{1}{2} \left( M \log(2\pi) + \log|R| + \nu_t^\top R^{-1} \nu_t \right) \right]\]

where $\nu_t = y_t - g(f(y_{t-1}, w_t))$ is the innovation (prediction error), $R$ is the observation noise covariance, and $M$ is the observation dimension.

The key difference from DirectIteration is that at each step the state is clamped to the observation: the prediction uses $f(y_{t-1}, \ldots)$ rather than $f(f(\ldots, u_0), \ldots)$.

AR(1) Example

using DifferenceEquations, LinearAlgebra, Random

rho_true = 0.8
sigma_e = 0.5
T = 200

# Simulate AR(1) data using the package
Random.seed!(42)
prob_sim = LinearStateSpaceProblem(
    fill(rho_true, 1, 1), fill(sigma_e, 1, 1), [0.0], (0, T))
sol_sim = solve(prob_sim)
y = sol_sim.u[2:end]  # observed states y_1, ..., y_T

# Compute conditional log-likelihood
prob = LinearStateSpaceProblem(
    fill(rho_true, 1, 1), nothing, [0.0], (0, T);
    observables = y,
    observables_noise = Diagonal([sigma_e^2]),
)
sol = solve(prob, ConditionalLikelihood())
sol.logpdf

-133.30350182762825

VAR(1) Example

The same approach works for multivariate models:

A = [0.8 0.1; -0.1 0.7]
B = [0.5 0.0; 0.0 0.5]
T_var = 100

# Simulate VAR(1) data
Random.seed!(123)
prob_sim_var = LinearStateSpaceProblem(A, B, zeros(2), (0, T_var))
sol_sim_var = solve(prob_sim_var)
y_var = sol_sim_var.u[2:end]

prob_var = LinearStateSpaceProblem(
    A, nothing, zeros(2), (0, T_var);
    observables = y_var,
    observables_noise = Diagonal([0.25, 0.25]),
)
sol_var = solve(prob_var, ConditionalLikelihood())
sol_var.logpdf

-148.67922110038018

Nonlinear Example with StateSpaceProblem

ConditionalLikelihood works with all problem types, including user-defined nonlinear callbacks via StateSpaceProblem.

Here we estimate a nonlinear AR(1): $x_{t+1} = \rho x_t + \alpha x_t^2 + e_t$. We first simulate data using a generic StateSpaceProblem, then compute the conditional likelihood.

rho_nl = 0.8
alpha_nl = 0.05
sigma_nl = 0.3
T_nl = 100

# Nonlinear transition (supports both mutable and immutable arrays)
function nl_transition!!(x_next, x, w, p, t)
    (; rho, alpha) = p
    val = rho * x[1] + alpha * x[1]^2
    if ismutable(x_next)
        x_next[1] = val + w[1]
        return x_next
    else
        return typeof(x)(val + w[1])
    end
end

p_nl = (; rho = rho_nl, alpha = alpha_nl)

# Simulate nonlinear data
Random.seed!(99)
prob_sim_nl = StateSpaceProblem(
    nl_transition!!, nothing, [0.0], (0, T_nl), p_nl;
    n_shocks = 1, n_obs = 0)
sol_sim_nl = solve(prob_sim_nl)
y_nl = sol_sim_nl.u[2:end]

# Conditional likelihood (no noise in prediction, noise only in obs)
function nl_transition_noiseless!!(x_next, x, w, p, t)
    (; rho, alpha) = p
    val = rho * x[1] + alpha * x[1]^2
    if ismutable(x_next)
        x_next[1] = val
        return x_next
    else
        return typeof(x)(val)
    end
end

prob_nl = StateSpaceProblem(
    nl_transition_noiseless!!, nothing, [0.0], (0, T_nl), p_nl;
    n_shocks = 0, n_obs = 0,
    observables = y_nl,
    observables_noise = Diagonal([sigma_nl^2]),
)
sol_nl = solve(prob_nl, ConditionalLikelihood())
sol_nl.logpdf

NaN

Maximum Likelihood Estimation

ConditionalLikelihood is fully differentiable with ForwardDiff.jl, making it straightforward to use with gradient-based optimization for MLE.

Use save_everystep=false when you only need the log-likelihood (not the full trajectory). This reduces allocations and speeds up ForwardDiff gradient computation — up to 7x faster with StaticArrays.

using ForwardDiff

# Negative log-likelihood as a function of rho
function neg_loglik(rho_vec)
    T_el = eltype(rho_vec)
    A_opt = fill(rho_vec[1], 1, 1)
    prob_opt = LinearStateSpaceProblem(
        A_opt, nothing, [zero(T_el)], (0, length(y));
        observables = y,
        observables_noise = Diagonal([T_el(sigma_e^2)]),
    )
    return -solve(prob_opt, ConditionalLikelihood(); save_everystep = false).logpdf
end

# Gradient at the true value
grad = ForwardDiff.gradient(neg_loglik, [rho_true])

1-element Vector{Float64}:
 6.106944176132063

The gradient is near zero at the true parameter value, confirming the MLE is correctly identified. For full optimization, use Optimization.jl with AutoForwardDiff() and an optimizer like LBFGS().

Using the Workspace

For repeated solves (e.g., inside an optimizer), use the init/solve! pattern to avoid repeated memory allocation:

ws = init(prob, ConditionalLikelihood())
sol_ws = solve!(ws)
sol_ws.logpdf

-133.30350182762825

With save_everystep=false, the workspace allocates only 2-element buffers:

ws_ep = init(prob, ConditionalLikelihood(); save_everystep = false)
sol_ep = solve!(ws_ep)
length(sol_ep.u)  # 2: [u_initial, u_final]

2