API Functions

diagnose(d::AbstractMacroData) -> DataDiagnostic

Scan data for NaN, Inf, constant columns, and very short series.

Examples

d = TimeSeriesData(randn(100, 3))
diag = diagnose(d)
diag.is_clean  # true if no issues

fix(d::TimeSeriesData; method=:listwise) -> TimeSeriesData

Fix data issues and return a clean copy.

Methods

• :listwise — drop rows with any NaN or Inf (default)
• :interpolate — linear interpolation for interior NaN, forward-fill edges
• :mean — replace NaN with column mean of finite values

Inf values are always replaced with NaN first (then handled by the chosen method). Constant columns are dropped with a warning.

Examples

d = TimeSeriesData([1.0 NaN; 2.0 3.0; 3.0 4.0])
d_clean = fix(d; method=:listwise)  # drops row 1

fix(d::PanelData; method=:listwise) -> PanelData

Fix data issues and return a clean copy.

Methods

• :listwise — drop rows with any NaN or Inf (default)
• :interpolate — linear interpolation within each group
• :mean — replace NaN with within-group column mean

Inf values are replaced with NaN first. Constant columns are dropped with a warning.

fix(d::CrossSectionData; method=:listwise) -> CrossSectionData

Fix data issues and return a clean copy.

Methods

• :listwise — drop rows with any NaN or Inf (default)
• :mean — replace NaN with column mean of finite values
• :interpolate — same as :mean (no time ordering in cross-section data)

Inf values are replaced with NaN first. Constant columns are dropped with a warning.

validate_for_model(d::AbstractMacroData, model_type::Symbol)

Check that data is compatible with the specified model type. Throws ArgumentError on mismatch.

Model types requiring multivariate data (n_vars ≥ 2)

:var, :vecm, :bvar, :factors, :dynamic_factors, :gdfm

Model types requiring univariate data (n_vars == 1)

:arima, :ar, :ma, :arma, :arch, :garch, :egarch, :gjr_garch, :sv, :hp_filter, :hamilton_filter, :beveridge_nelson, :baxter_king, :boosted_hp, :adf, :kpss, :pp, :za, :ngperron

Model types accepting any dimensionality

:lp, :lp_iv, :smooth_lp, :state_lp, :propensity_lp, :gmm

Examples

d = TimeSeriesData(randn(100, 3))
validate_for_model(d, :var)    # OK
validate_for_model(d, :arima)  # throws ArgumentError

apply_tcode(y::AbstractVector{<:Real}, tcode::Int) -> Vector{Float64}

Apply FRED transformation code to a univariate series.

Codes 4–7 require strictly positive data. The output vector is shorter than the input for difference-based codes:

• tcode 1: same length
• tcode 2, 4, 5: length T-1
• tcode 3, 6, 7: length T-2

Examples

y = [100.0, 102.0, 105.0, 103.0, 108.0]
apply_tcode(y, 5)  # log first differences (approx growth rates)

apply_tcode(d::TimeSeriesData, tcodes::Vector{Int}) -> TimeSeriesData

Apply per-variable FRED transformation codes. Rows are trimmed consistently (to the shortest transformed series).

Examples

d = TimeSeriesData(rand(200, 3) .+ 1; varnames=["GDP","CPI","FFR"])
d2 = apply_tcode(d, [5, 5, 1])  # log-diff GDP and CPI, leave FFR in levels

apply_tcode(d::TimeSeriesData, tcode::Int) -> TimeSeriesData

Apply the same FRED transformation code to all variables.

apply_tcode(d::TimeSeriesData) -> TimeSeriesData

Convenience method: apply the transformation codes stored in d.tcode. Equivalent to apply_tcode(d, d.tcode).

Examples

d = load_example(:fred_md)  # tcode already set
d_transformed = apply_tcode(d)

apply_tcode(d::PanelData, tcodes::Vector{Int}) -> PanelData

Apply per-variable FRED transformation codes to a PanelData container, processing each group independently. Each group is trimmed to its own common length after transformation, then reassembled.

Examples

pd = xtset(df, :id, :t)
pd2 = apply_tcode(pd, [5, 5, 1])  # log-diff first two vars, leave third in levels

apply_tcode(d::PanelData, tcode::Int) -> PanelData

Apply the same FRED transformation code to all variables in a PanelData container.

inverse_tcode(y::AbstractVector{<:Real}, tcode::Int;
              x_prev::Union{AbstractVector{<:Real},Nothing}=nothing) -> Vector{Float64}

Undo a FRED transformation to recover the original series.

For difference-based codes (2, 3, 5, 6, 7), x_prev must supply the initial values needed for reconstruction:

• tcode 2: x_prev = [x₀] (1 value)
• tcode 3: x_prev = [x₀, x₁] (2 values, original levels)
• tcode 5: x_prev = [x₀] (1 value, original level)
• tcode 6: x_prev = [x₀, x₁] (2 values, original levels)
• tcode 7: x_prev = [x₀, x₁] (2 values, original levels)

Examples

y = [100.0, 102.0, 105.0]
yd = apply_tcode(y, 2)                        # [2.0, 3.0]
inverse_tcode(yd, 2; x_prev=[100.0])          # [102.0, 105.0]

apply_filter(d::TimeSeriesData, specs::AbstractVector; component=:cycle, kwargs...)

Apply per-variable filter specifications to a TimeSeriesData container.

Each element of specs can be:

• Symbol — filter name (:hp, :hamilton, :bn, :bk, :boosted_hp)
• AbstractFilterResult — a pre-computed filter result
• Tuple{filter, Symbol} — filter with per-variable component override (e.g., (:hp, :trend))
• nothing — pass-through (no filtering for this variable)

The output is trimmed to the intersection of valid ranges across all variables.

Arguments

• d::TimeSeriesData — input data
• specs::AbstractVector — per-variable filter specifications (length must equal n_vars)

Keyword Arguments

• component::Symbol=:cycle — default component to extract (:cycle or :trend)
• Additional kwargs are forwarded to filter functions

Returns

A new TimeSeriesData with filtered data, trimmed to the common valid range.

Examples

d = TimeSeriesData(cumsum(randn(200, 3), dims=1); varnames=["GDP","CPI","FFR"])

# Per-variable: HP cycle for GDP, Hamilton cycle for CPI, pass-through FFR
d2 = apply_filter(d, [:hp, :hamilton, nothing])

# Per-variable component overrides via tuples
d3 = apply_filter(d, [(:hp, :trend), (:hamilton, :cycle), nothing])

apply_filter(d::TimeSeriesData, spec::Union{Symbol, AbstractFilterResult};
             component=:cycle, vars=nothing, kwargs...)

Apply a single filter to all variables (or a subset specified by vars).

Variables not in vars are passed through unchanged.

Arguments

• d::TimeSeriesData — input data
• spec — filter symbol (:hp, :hamilton, :bn, :bk, :boosted_hp) or pre-computed result

Keyword Arguments

• component::Symbol=:cycle — component to extract (:cycle or :trend)
• vars::Union{Nothing, Vector{String}, Vector{Int}} — variables to filter (default: all)
• Additional kwargs are forwarded to filter functions

Examples

d = TimeSeriesData(cumsum(randn(200, 3), dims=1); varnames=["GDP","CPI","FFR"])

# HP cycle for all variables
d_hp = apply_filter(d, :hp)

# HP trend for selected variables only
d_sel = apply_filter(d, :hp; vars=["GDP", "CPI"], component=:trend)

apply_filter(d::PanelData, spec; component=:cycle, vars=nothing, kwargs...)

Apply filters to a PanelData container group-by-group.

Each group is extracted via group_data, filtered, and the results are reassembled into a new PanelData. Each group is trimmed independently to its own common valid range (groups may have different resulting lengths if unbalanced).

Arguments

• d::PanelData — input panel data
• spec — filter specification: a single Symbol, AbstractFilterResult, or AbstractVector of per-variable specs (same formats as TimeSeriesData)

Keyword Arguments

• component::Symbol=:cycle — component to extract
• vars — variables to filter (others pass-through)
• Additional kwargs forwarded to filter functions

Examples

pd = xtset(df, :id, :t)
pd_hp = apply_filter(pd, :hp; component=:cycle)
pd_sel = apply_filter(pd, :hp; vars=["GDP"], component=:trend)

describe_data(d::AbstractMacroData) -> DataSummary

Compute per-variable summary statistics (N, Mean, Std, Min, P25, Median, P75, Max, Skewness, Kurtosis). NaN values are excluded from all computations.

For PanelData, also prints panel dimensions via panel_summary.

Examples

d = TimeSeriesData(randn(200, 3); varnames=["GDP","CPI","FFR"], frequency=Quarterly)
describe_data(d)

xtset(df::DataFrame, group_col::Symbol, time_col::Symbol;
      varnames=nothing, frequency=Other, tcode=nothing,
      cohort=nothing) -> PanelData{Float64}

Construct a PanelData container from a DataFrame, analogous to Stata's xtset.

Extracts all numeric columns (excluding group_col, time_col, and cohort), sorts by (group, time), validates no duplicate (group, time) pairs, and detects whether the panel is balanced.

Arguments

• df::DataFrame — input data
• group_col::Symbol — column name identifying groups (entities)
• time_col::Symbol — column name identifying time periods

Keyword Arguments

• varnames::Union{Vector{String},Nothing} — override variable names (default: column names)
• frequency::Frequency — data frequency (default: Other)
• tcode::Union{Vector{Int},Nothing} — transformation codes per variable
• cohort::Union{Symbol,Nothing} — column name identifying cohort membership (default: nothing)

Examples

using DataFrames
df = DataFrame(id=repeat(1:3, inner=50), t=repeat(1:50, 3),
               x=randn(150), y=randn(150))
pd = xtset(df, :id, :t)

isbalanced(d::PanelData) -> Bool

Return true if all groups have the same number of observations.

groups(d::PanelData) -> Vector{String}

Return the group names.

ngroups(d::PanelData) -> Int

Return the number of groups.

group_data(d::PanelData, g) -> TimeSeriesData

Extract data for a single group as a TimeSeriesData container.

g can be an integer group index or a string group name.

Examples

pd = xtset(df, :id, :t)
g1 = group_data(pd, 1)       # by index
g1 = group_data(pd, "1")     # by name

panel_summary(d::PanelData)

Display a summary table of the panel structure: number of groups, observations per group (min, mean, max), and balance status.

Examples

pd = xtset(df, :id, :t)
panel_summary(pd)

StatsAPI.nobs(d::AbstractMacroData)

Return the number of observations.

nvars(d::AbstractMacroData)

Return the number of variables.

varnames(d::AbstractMacroData)

Return variable names.

frequency(d::TimeSeriesData)
frequency(d::PanelData)

Return the data frequency.

time_index(d::TimeSeriesData)

Return the integer time index vector.

obs_id(d::CrossSectionData)

Return the observation identifier vector.

desc(d::AbstractMacroData) -> String

Return the dataset description. Returns "" if no description has been set.

Examples

d = TimeSeriesData(randn(100, 3); desc="US macro quarterly data")
desc(d)  # "US macro quarterly data"

vardesc(d::AbstractMacroData) -> Dict{String,String}

Return the dictionary of per-variable descriptions.

Examples

d = TimeSeriesData(randn(100, 2); varnames=["GDP","CPI"],
    vardesc=Dict("GDP" => "Real GDP growth", "CPI" => "Consumer prices"))
vardesc(d)  # Dict("GDP" => "Real GDP growth", "CPI" => "Consumer prices")

vardesc(d::AbstractMacroData, name::String) -> String

Return the description for variable name. Throws ArgumentError if no description exists for that variable.

Examples

vardesc(d, "GDP")  # "Real GDP growth"

rename_vars!(d::AbstractMacroData, old => new)
rename_vars!(d::AbstractMacroData, names::Vector{String})

Rename variables in a data container. With a Pair, renames a single variable. With a Vector{String}, replaces all variable names. Also updates vardesc keys.

set_time_index!(d::TimeSeriesData, idx::Vector{Int})

Set the time index for a TimeSeriesData container.

set_obs_id!(d::CrossSectionData, ids::Vector{Int})

Set observation identifiers for a CrossSectionData container.

set_desc!(d::AbstractMacroData, text::String)

Set the dataset description.

Examples

d = TimeSeriesData(randn(100, 3))
set_desc!(d, "US macroeconomic quarterly data 1959-2024")
desc(d)  # "US macroeconomic quarterly data 1959-2024"

set_vardesc!(d::AbstractMacroData, name::String, text::String)

Set the description for a single variable. The variable must exist in the data.

Examples

d = TimeSeriesData(randn(100, 2); varnames=["GDP", "CPI"])
set_vardesc!(d, "GDP", "Real Gross Domestic Product, seasonally adjusted")
set_vardesc!(d, "CPI", "Consumer Price Index for All Urban Consumers")
vardesc(d, "GDP")  # "Real Gross Domestic Product, seasonally adjusted"

set_vardesc!(d::AbstractMacroData, descriptions::Dict{String,String})

Set descriptions for multiple variables at once. All keys must be valid variable names.

Examples

set_vardesc!(d, Dict("GDP" => "Real GDP", "CPI" => "Consumer prices"))

dates(d::TimeSeriesData) -> Vector{String}

Return the date labels vector. Returns empty vector if no dates are set.

set_dates!(d::TimeSeriesData, dates::Vector{String})

Set the date labels for a TimeSeriesData container. Length must match T_obs, or pass an empty vector to clear dates.

Examples

d = TimeSeriesData(randn(4, 2))
set_dates!(d, ["2020Q1", "2020Q2", "2020Q3", "2020Q4"])
dates(d)  # ["2020Q1", "2020Q2", "2020Q3", "2020Q4"]

to_matrix(d::TimeSeriesData) -> Matrix

Return the raw data matrix from a TimeSeriesData container.

to_matrix(d::PanelData) -> Matrix

Return the raw stacked data matrix from a PanelData container.

to_matrix(d::CrossSectionData) -> Matrix

Return the raw data matrix from a CrossSectionData container.

to_vector(d::TimeSeriesData) -> Vector

Return the data as a vector (requires exactly 1 variable).

to_vector(d::TimeSeriesData, var::Int) -> Vector

Return a single column by index.

to_vector(d::TimeSeriesData, var::String) -> Vector

Return a single column by name.

load_example(name::Symbol) -> AbstractMacroData

Load a built-in example dataset.

Available Datasets

• :fred_md — FRED-MD Monthly Database, January 2026 vintage (126 variables × 804 months) → TimeSeriesData
• :fred_qd — FRED-QD Quarterly Database, January 2026 vintage (245 variables × 268 quarters) → TimeSeriesData
• :pwt — Penn World Table 10.01, 38 OECD countries (42 variables × 74 years, 1950–2023) → PanelData
• :ddcg — Acemoglu et al. (2019) Democracy-GDP Panel (2 variables × 9,384 obs, 184 countries × 51 years) → PanelData
• :mpdta — Callaway & Sant'Anna (2021) Minimum Wage Panel (3 variables × 2,500 obs, 500 counties × 5 years) → PanelData

For time series datasets, the returned TimeSeriesData includes variable names, transformation codes, frequency, per-variable descriptions (via vardesc), dataset description (via desc), and bibliographic references (via refs).

For panel datasets, the returned PanelData includes country identifiers as groups, year identifiers as time index, variable descriptions, and references.

Examples

# Load FRED-MD
md = load_example(:fred_md)
nobs(md)       # 804
nvars(md)      # 126
desc(md)       # "FRED-MD Monthly Database, January 2026 Vintage (McCracken & Ng 2016)"
vardesc(md, "INDPRO")  # "IP Index"
refs(md)       # McCracken & Ng (2016)

# Apply recommended transformations
md_transformed = apply_tcode(md, md.tcode)

# Load FRED-QD
qd = load_example(:fred_qd)

# Load Penn World Table (panel data)
pwt = load_example(:pwt)
nobs(pwt)         # 2812 (38 countries × 74 years)
nvars(pwt)        # 42
ngroups(pwt)      # 38
groups(pwt)       # ["AUS", "AUT", ..., "USA"]
isbalanced(pwt)   # true
g = group_data(pwt, "USA")  # extract single country as TimeSeriesData
refs(pwt)         # Feenstra, Inklaar & Timmer (2015)

hp_filter(y::AbstractVector; lambda=1600.0) -> HPFilterResult

Apply the Hodrick-Prescott filter to decompose time series y into trend and cycle.

Solves the optimization problem:

\[\min_\tau \sum_{t=1}^T (y_t - \tau_t)^2 + \lambda \sum_{t=2}^{T-1} (\tau_{t+1} - 2\tau_t + \tau_{t-1})^2\]

Implementation uses a sparse pentadiagonal Cholesky factorization for O(T) cost.

Arguments

• y::AbstractVector: Time series data (length ≥ 3)

Keywords

• lambda::Real=1600.0: Smoothing parameter. Common values: 6.25 (annual), 1600 (quarterly, default), 129600 (monthly)

Returns

• HPFilterResult{T} with fields trend, cycle, lambda, T_obs

Examples

y = cumsum(randn(200))
result = hp_filter(y)
result = hp_filter(y; lambda=6.25)  # annual data

References

• Hodrick, R. J., & Prescott, E. C. (1997). JMCB 29(1): 1–16.

hamilton_filter(y::AbstractVector; h=8, p=4) -> HamiltonFilterResult

Apply the Hamilton (2018) regression filter for trend-cycle decomposition.

Regresses $y_{t+h}$ on $[1, y_t, y_{t-1}, \ldots, y_{t-p+1}]$ by OLS. The residuals form the cyclical component and the fitted values form the trend.

The filter loses h + p - 1 observations at the start of the sample.

Arguments

• y::AbstractVector: Time series data

Keywords

• h::Int=8: Forecast horizon (default 8 for quarterly data = 2 years)
• p::Int=4: Number of lags in the regression (default 4 for quarterly data)

Returns

• HamiltonFilterResult{T} with fields trend, cycle, beta, h, p, T_obs, valid_range

Examples

y = cumsum(randn(200))
result = hamilton_filter(y)                 # quarterly defaults
result = hamilton_filter(y; h=24, p=12)     # monthly data (2-year horizon)

References

• Hamilton, J. D. (2018). REStat 100(5): 831–843.

beveridge_nelson(y::AbstractVector; method=:arima, p=:auto, q=:auto, max_terms=500, cycle_order=2) -> BeveridgeNelsonResult

Compute the Beveridge-Nelson decomposition of time series y.

Assumes y is I(1) and decomposes it into:

\[y_t = \tau_t + c_t\]

where $\tau_t$ is a random walk with drift (permanent component) and $c_t$ is a stationary transitory component.

Methods

• :arima (default): Classic BN via ARIMA representation (Beveridge & Nelson, 1981)
• :statespace: Correlated UC model via MLE + Kalman smoother (Morley, Nelson & Zivot, 2003)

Arguments

• y::AbstractVector: Time series data (assumed I(1), length ≥ 10)

Keywords

• method::Symbol=:arima: Decomposition method
• p: AR order for ARMA model of Δy (method=:arima). :auto uses auto_arima
• q: MA order for ARMA model of Δy (method=:arima). :auto uses auto_arima
• max_terms::Int=500: Maximum ψ-weights for MA(∞) truncation (method=:arima)
• cycle_order::Int=2: AR order for cyclical component (method=:statespace, 1 or 2)

Returns

• BeveridgeNelsonResult{T} with fields permanent, transitory, drift, long_run_multiplier, arima_order, T_obs

Examples

# Classic ARIMA-based BN decomposition
y = cumsum(randn(200)) + 0.3 * sin.(2π * (1:200) / 20)
result = beveridge_nelson(y)
result = beveridge_nelson(y; p=2, q=1)  # manual ARMA order

# Morley (2002) correlated UC model
result = beveridge_nelson(y; method=:statespace)
result = beveridge_nelson(y; method=:statespace, cycle_order=1)

References

• Beveridge, S., & Nelson, C. R. (1981). JME 7(2): 151–174.
• Morley, J. C., Nelson, C. R., & Zivot, E. (2003). REStat 85(2): 235–243.

baxter_king(y::AbstractVector; pl=6, pu=32, K=12) -> BaxterKingResult

Apply the Baxter-King band-pass filter to isolate business cycle frequencies.

Computes symmetric moving average weights that approximate the ideal band-pass filter for periods between pl and pu. The filter is constrained to sum to zero (removing stochastic trends).

Loses K observations at each end of the sample (2K total).

Arguments

• y::AbstractVector: Time series data (length > 2K)

Keywords

• pl::Int=6: Minimum period of oscillation to pass (quarterly: 6 = 1.5 years)
• pu::Int=32: Maximum period of oscillation to pass (quarterly: 32 = 8 years)
• K::Int=12: Truncation length (number of leads/lags in the moving average)

Returns

• BaxterKingResult{T} with fields cycle, trend, weights, pl, pu, K, T_obs, valid_range

Examples

y = cumsum(randn(200))
result = baxter_king(y)                    # quarterly defaults (6-32 quarters)
result = baxter_king(y; pl=2, pu=8, K=6)  # annual data (2-8 years)

References

• Baxter, M., & King, R. G. (1999). REStat 81(4): 575–593.

boosted_hp(y::AbstractVector; lambda=1600.0, stopping=:BIC, max_iter=100, sig_p=0.05) -> BoostedHPResult

Apply the boosted HP filter (Phillips & Shi 2021) for improved trend estimation.

Iteratively re-filters the cyclical component of the standard HP filter. At each iteration, the cycle is decomposed again and the newly estimated "trend of the cycle" is added back to the overall trend.

Stopping criteria

• :ADF — Stop when the ADF test rejects the null of a unit root in the cycle at significance level sig_p (recommended for detecting stationarity)
• :BIC — Stop when the Phillips-Shi information criterion increases (balances variance reduction against effective degrees of freedom)
• :fixed — Run all max_iter iterations

Arguments

• y::AbstractVector: Time series data (length ≥ 3)

Keywords

• lambda::Real=1600.0: HP smoothing parameter
• stopping::Symbol=:BIC: Stopping criterion (:ADF, :BIC, or :fixed)
• max_iter::Int=100: Maximum number of boosting iterations
• sig_p::Real=0.05: Significance level for ADF stopping criterion

Returns

• BoostedHPResult{T} with fields trend, cycle, lambda, iterations, stopping, bic_path, adf_pvalues, T_obs

Examples

y = cumsum(randn(200))
result = boosted_hp(y)                              # BIC stopping (default)
result = boosted_hp(y; stopping=:ADF, sig_p=0.05)   # ADF stopping
result = boosted_hp(y; stopping=:fixed, max_iter=5)  # fixed iterations

References

• Phillips, P. C. B., & Shi, Z. (2021). IER 62(2): 521–570.
• Mei, Z., Phillips, P. C. B., & Shi, Z. (2024). JAE 39(7): 1260–1281.

trend(result::AbstractFilterResult) -> Vector

Return the trend component from a filter result. For BeveridgeNelsonResult, returns the permanent component.

cycle(result::AbstractFilterResult) -> Vector

Return the cyclical component from a filter result. For BeveridgeNelsonResult, returns the transitory component.

estimate_ar(y, p; method=:ols, include_intercept=true) -> ARModel

Estimate AR(p) model: yₜ = c + φ₁yₜ₋₁ + ... + φₚyₜ₋ₚ + εₜ

Arguments

• y: Time series vector
• p: AR order (must be ≥ 1)
• method: Estimation method (:ols or :mle)
• include_intercept: Whether to include constant term

Returns

ARModel with estimated coefficients and diagnostics.

Example

y = randn(200)
model = estimate_ar(y, 2)
println(model.phi)  # AR coefficients

estimate_ma(y, q; method=:css_mle, include_intercept=true, max_iter=500) -> MAModel

Estimate MA(q) model: yₜ = c + εₜ + θ₁εₜ₋₁ + ... + θqεₜ₋q

Arguments

• y: Time series vector
• q: MA order (must be ≥ 1)
• method: Estimation method (:css, :mle, or :css_mle)
• include_intercept: Whether to include constant term
• max_iter: Maximum optimization iterations

Returns

MAModel with estimated coefficients and diagnostics.

Example

y = randn(200)
model = estimate_ma(y, 1)
println(model.theta)  # MA coefficient

estimate_arma(y, p, q; method=:css_mle, include_intercept=true, max_iter=500) -> ARMAModel

Estimate ARMA(p,q) model: yₜ = c + φ₁yₜ₋₁ + ... + φₚyₜ₋ₚ + εₜ + θ₁εₜ₋₁ + ... + θqεₜ₋q

Arguments

• y: Time series vector
• p: AR order
• q: MA order
• method: Estimation method (:css, :mle, or :css_mle)
• include_intercept: Whether to include constant term
• max_iter: Maximum optimization iterations

Returns

ARMAModel with estimated coefficients and diagnostics.

Example

y = randn(200)
model = estimate_arma(y, 1, 1)
println("AR: ", model.phi, " MA: ", model.theta)

estimate_arima(y, p, d, q; method=:css_mle, include_intercept=true, max_iter=500) -> ARIMAModel

Estimate ARIMA(p,d,q) model by differencing d times and fitting ARMA(p,q).

Arguments

• y: Time series vector
• p: AR order
• d: Integration order (number of differences)
• q: MA order
• method: Estimation method (:css, :mle, or :css_mle)
• include_intercept: Whether to include constant term (on differenced series)
• max_iter: Maximum optimization iterations

Returns

ARIMAModel with estimated coefficients and diagnostics.

Example

y = cumsum(randn(200))  # Random walk
model = estimate_arima(y, 1, 1, 0)  # ARIMA(1,1,0)
println(model.phi)

_compute_psi_weights(phi, theta, h) -> Vector{T}

Compute ψ-weights for the MA(∞) representation of an ARMA process.

The ARMA(p,q) process can be written as: yₜ = μ + Σⱼ₌₀^∞ ψⱼ εₜ₋ⱼ

where ψ₀ = 1 and ψⱼ follows the recursion: ψⱼ = φ₁ψⱼ₋₁ + ... + φₚψⱼ₋ₚ + θⱼ

Returns ψ₁, ψ₂, ..., ψₕ.

_confidence_band(forecasts, se, conf_level)

Compute symmetric confidence interval bounds from forecasts, standard errors, and a confidence level.

_forecast_arma(y, resid, c, phi, theta, sigma2, h, conf_level) -> ARIMAForecast

Unified point forecast + CI computation for any ARMA(p,q) model. AR models pass theta=T[], MA models pass phi=T[].

_forecast_variance(sigma2, psi, h) -> Vector{T}

Compute h-step ahead forecast variance.

Var(eₜ₊ₕ) = σ² (1 + ψ₁² + ψ₂² + ... + ψₕ₋₁²)

_integrate_forecasts(y, fc_diff, d) -> Vector{T}

Integrate d-differenced forecasts back to original scale.

_integrate_se(se_diff, d) -> Vector{T}

Approximate standard errors after integration.

For d-fold integration, the variance grows roughly as h^d. This is a conservative approximation.

forecast(model::ARIMAModel, h; conf_level=0.95) -> ARIMAForecast

Compute h-step ahead forecasts with confidence intervals for ARIMA model. Forecasts are computed on the differenced series and then integrated back to the original scale.

forecast(model::ARMAModel, h; conf_level=0.95) -> ARIMAForecast

Compute h-step ahead forecasts with confidence intervals for ARMA model.

forecast(model::ARModel, h; conf_level=0.95) -> ARIMAForecast

Compute h-step ahead forecasts with confidence intervals for AR model.

forecast(model::MAModel, h; conf_level=0.95) -> ARIMAForecast

Compute h-step ahead forecasts with confidence intervals for MA model.

predict(model::AbstractARIMAModel, h::Int) -> Vector{T}

Return h-step ahead point forecasts (without confidence intervals).

select_arima_order(y, max_p, max_q; criterion=:bic, d=0, method=:css_mle, include_intercept=true)

Automatically select ARMA/ARIMA order via grid search over information criteria.

Searches over p ∈ 0:maxp and q ∈ 0:maxq, fits each model, and selects the order that minimizes the specified information criterion.

Arguments

• y: Time series vector
• max_p: Maximum AR order to consider
• max_q: Maximum MA order to consider
• criterion: Selection criterion (:aic or :bic, default :bic)
• d: Integration order for ARIMA (default 0 = ARMA)
• method: Estimation method (:css, :mle, or :css_mle)
• include_intercept: Whether to include constant term

Returns

ARIMAOrderSelection with best orders, IC matrices, and fitted models.

Example

y = randn(200)
result = select_arima_order(y, 3, 3; criterion=:bic)
println("Best order: p=$(result.best_p_bic), q=$(result.best_q_bic)")
best_model = result.best_model_bic

auto_arima(y; max_p=5, max_q=5, max_d=2, criterion=:bic, method=:css_mle, stepwise=true)

Automatically select and fit the best ARIMA model.

Performs order selection over p, d, and q using the specified criterion. For integration order d, uses unit root test heuristics.

Arguments

• y: Time series vector
• max_p: Maximum AR order (default 5)
• max_q: Maximum MA order (default 5)
• max_d: Maximum integration order (default 2)
• criterion: Selection criterion (:aic or :bic)
• method: Estimation method
• stepwise: If true (default), use Hyndman-Khandakar (2008) stepwise search. If false, use exhaustive grid search over all (p,q) combinations.

Returns

Best fitted ARIMAModel or ARMAModel.

Example

y = cumsum(randn(200))
model = auto_arima(y)
println(model)

References

• Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: the forecast package for R. Journal of Statistical Software 27(3): 1–22.

ic_table(result::ARIMAOrderSelection; criterion=:bic)

Return a formatted table of IC values for printing.

Return AR order p.

Return MA order q.

Return integration order d.

Standard errors for AR model coefficients.

Standard errors for MA model coefficients.

Standard errors for ARMA model coefficients.

Standard errors for ARIMA model coefficients.

estimate_reg(y, X; cov_type=:hc1, weights=nothing, varnames=nothing, clusters=nothing) -> RegModel{T}

Estimate a linear regression model via OLS or WLS.

OLS: beta = (X'X)^{-1} X'y WLS: beta = (X'WX)^{-1} X'Wy where W = diag(weights)

Arguments

• y::AbstractVector{T} — dependent variable (n x 1)
• X::AbstractMatrix{T} — regressor matrix (n x k), should include a constant column
• cov_type::Symbol — covariance estimator: :ols, :hc0, :hc1 (default), :hc2, :hc3, :cluster
• weights::Union{Nothing,AbstractVector} — WLS weights (nothing for OLS)
• varnames::Union{Nothing,Vector{String}} — coefficient names (auto-generated if nothing)
• clusters::Union{Nothing,AbstractVector} — cluster assignments (required for :cluster)

Returns

RegModel{T} with estimated coefficients, robust covariance matrix, fit statistics.

Examples

using MacroEconometricModels
n = 200
X = hcat(ones(n), randn(n, 2))
beta_true = [1.0, 2.0, -0.5]
y = X * beta_true + 0.5 * randn(n)
m = estimate_reg(y, X)
report(m)

References

• Greene, W. H. (2018). Econometric Analysis. 8th ed. Pearson.
• White, H. (1980). Econometrica 48(4), 817-838.

estimate_iv(y, X, Z; endogenous, cov_type=:hc1, varnames=nothing) -> RegModel{T}

Estimate a linear regression model via instrumental variables / two-stage least squares (IV/2SLS).

Algorithm

1. Project regressors onto instrument space: Xhat = PZ * X where P_Z = Z(Z'Z)^{-1}Z'
2. Second stage: beta = (Xhat'X)^{-1} Xhat'y
3. Residuals from ORIGINAL X: e = y - Xbeta (not X_hatbeta)
4. Covariance uses X_hat for the bread and original residuals for the meat

Arguments

• y::AbstractVector{T} — dependent variable (n x 1)
• X::AbstractMatrix{T} — regressor matrix (n x k), includes exogenous regressors and intercept
• Z::AbstractMatrix{T} — instrument matrix (n x m), includes exogenous regressors and excluded instruments
• endogenous::Vector{Int} — column indices of endogenous regressors in X
• cov_type::Symbol — covariance estimator: :ols, :hc0, :hc1 (default), :hc2, :hc3
• varnames::Union{Nothing,Vector{String}} — coefficient names (auto-generated if nothing)

Returns

RegModel{T} with method=:iv, including first-stage F-statistic and Sargan test.

Diagnostics

• first_stage_f: minimum first-stage F across endogenous variables (>10 suggests strong instruments)
• sargan_stat/sargan_pval: overidentification test (nothing if exactly identified)

Examples

using MacroEconometricModels
n = 500
z1, z2 = randn(n), randn(n)
x_endog = 0.5 * z1 + 0.3 * z2 + randn(n)
u = randn(n)
x_endog .+= 0.5 * u  # endogeneity
y = 1.0 .+ 2.0 * x_endog + u
X = hcat(ones(n), x_endog)
Z = hcat(ones(n), z1, z2)
m = estimate_iv(y, X, Z; endogenous=[2])
report(m)

References

• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press, ch. 5.
• Stock, J. H. & Yogo, M. (2005). Identification and Inference for Econometric Models. Cambridge University Press, ch. 5.

estimate_logit(y, X; cov_type=:ols, varnames=nothing, clusters=nothing, maxiter=100, tol=1e-8) -> LogitModel{T}

Estimate a binary logistic regression model via maximum likelihood (IRLS/Fisher scoring).

Algorithm

Uses iteratively reweighted least squares with the canonical logit link mu = 1/(1+exp(-X*beta)). Converges when the change in log-likelihood is below tol * (|loglik| + 1).

Arguments

• y::AbstractVector{T} — binary dependent variable (0/1), n x 1
• X::AbstractMatrix{T} — regressor matrix (n x k), should include a constant column
• cov_type::Symbol — covariance estimator: :ols (information matrix), :hc0, :hc1, :hc2, :hc3, :cluster
• varnames::Union{Nothing,Vector{String}} — coefficient names (auto-generated if nothing)
• clusters::Union{Nothing,AbstractVector} — cluster assignments (required for :cluster)
• maxiter::Int — maximum IRLS iterations (default 100)
• tol — convergence tolerance (default 1e-8)

Returns

LogitModel{T} with estimated coefficients, McFadden pseudo-R-squared, deviance residuals.

Examples

using MacroEconometricModels, Random
rng = MersenneTwister(42)
n = 500
X = hcat(ones(n), randn(rng, n, 2))
beta_true = [0.0, 1.5, -1.0]
p = 1 ./ (1 .+ exp.(-X * beta_true))
y = Float64.(rand(rng, n) .< p)
m = estimate_logit(y, X; varnames=["const", "x1", "x2"])
report(m)

References

• McCullagh, P. & Nelder, J. A. (1989). Generalized Linear Models. Chapman & Hall.
• Agresti, A. (2002). Categorical Data Analysis. 2nd ed. Wiley.

estimate_probit(y, X; cov_type=:ols, varnames=nothing, clusters=nothing, maxiter=100, tol=1e-8) -> ProbitModel{T}

Estimate a binary probit regression model via maximum likelihood (IRLS/Fisher scoring).

Algorithm

Uses iteratively reweighted least squares with the probit link mu = Phi(X*beta), where Phi is the standard normal CDF. Converges when the change in log-likelihood is below tol * (|loglik| + 1).

Arguments

• y::AbstractVector{T} — binary dependent variable (0/1), n x 1
• X::AbstractMatrix{T} — regressor matrix (n x k), should include a constant column
• cov_type::Symbol — covariance estimator: :ols (information matrix), :hc0, :hc1, :hc2, :hc3, :cluster
• varnames::Union{Nothing,Vector{String}} — coefficient names (auto-generated if nothing)
• clusters::Union{Nothing,AbstractVector} — cluster assignments (required for :cluster)
• maxiter::Int — maximum IRLS iterations (default 100)
• tol — convergence tolerance (default 1e-8)

Returns

ProbitModel{T} with estimated coefficients, McFadden pseudo-R-squared, deviance residuals.

Examples

using MacroEconometricModels, Random
rng = MersenneTwister(42)
n = 500
X = hcat(ones(n), randn(rng, n, 2))
beta_true = [0.0, 1.0, -0.8]
using Distributions
p = cdf.(Normal(), X * beta_true)
y = Float64.(rand(rng, n) .< p)
m = estimate_probit(y, X; varnames=["const", "x1", "x2"])
report(m)

References

• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press, ch. 15.
• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 14.

marginal_effects(m::LogitModel{T}; type=:ame, at=nothing, conf_level=0.95) -> MarginalEffects{T}

Compute marginal effects from a logistic regression model with delta-method SEs.

Types

• :ame — Average Marginal Effects: (1/N) sumi f(Xi beta) * beta_j
• :mem — Marginal Effects at the Mean: f(Xbar beta) * betaj
• :mer — Marginal Effects at Representative values: f(x0 beta) * betaj

where f(eta) = logistic PDF = p(1-p) and p = 1/(1+exp(-eta)).

Arguments

• m::LogitModel{T} — estimated logit model
• type::Symbol — :ame (default), :mem, or :mer
• at::Union{Nothing,Dict} — for :mer, Dict mapping column index => value
• conf_level::Real — confidence level (default 0.95)

Returns

MarginalEffects{T} with effects, delta-method SEs, z-stats, p-values, and CIs.

Examples

m = estimate_logit(y, X)
me_ame = marginal_effects(m)                          # AME
me_mem = marginal_effects(m; type=:mem)                # MEM
me_mer = marginal_effects(m; type=:mer, at=Dict(2=>0.0))  # MER at x2=0

References

• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 14.
• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press.

marginal_effects(m::ProbitModel{T}; type=:ame, at=nothing, conf_level=0.95) -> MarginalEffects{T}

Compute marginal effects from a probit regression model with delta-method SEs.

Types

• :ame — Average Marginal Effects: (1/N) sumi phi(Xi beta) * beta_j
• :mem — Marginal Effects at the Mean: phi(Xbar beta) * betaj
• :mer — Marginal Effects at Representative values: phi(x0 beta) * betaj

where phi(eta) = standard normal PDF.

Arguments

• m::ProbitModel{T} — estimated probit model
• type::Symbol — :ame (default), :mem, or :mer
• at::Union{Nothing,Dict} — for :mer, Dict mapping column index => value
• conf_level::Real — confidence level (default 0.95)

Returns

MarginalEffects{T} with effects, delta-method SEs, z-stats, p-values, and CIs.

Examples

m = estimate_probit(y, X)
me_ame = marginal_effects(m)                          # AME
me_mem = marginal_effects(m; type=:mem)                # MEM
me_mer = marginal_effects(m; type=:mer, at=Dict(2=>0.0))  # MER at x2=0

References

• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 14.

marginal_effects(m::OrderedLogitModel{T}) -> NamedTuple

Compute average marginal effects (AME) for an ordered logit model.

Returns a K x J matrix of AMEs where element (k,j) is the average marginal effect of variable k on P(y = j).

Formula: dP(y=j)/dxk = [f(alpha{j-1} - x'beta) - f(alphaj - x'beta)] * betak where f is the logistic PDF, alpha0 = -Inf, alphaJ = +Inf.

Key property: AMEs sum to zero across categories for each variable.

Returns

Named tuple with fields:

• effects::Matrix{T} – K x J matrix of AMEs
• varnames::Vector{String} – variable names
• categories::Vector – category labels

References

• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 15.
• Greene, W. H. (2012). Econometric Analysis. 7th ed. Prentice Hall.

marginal_effects(m::OrderedProbitModel{T}) -> NamedTuple

Compute average marginal effects (AME) for an ordered probit model.

Same structure as marginal_effects(::OrderedLogitModel) but using the standard normal PDF.

References

• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 15.
• Greene, W. H. (2012). Econometric Analysis. 7th ed. Prentice Hall.

marginal_effects(m::MultinomialLogitModel{T}) -> NamedTuple

Compute average marginal effects (AME) for a multinomial logit model.

Returns a K x J matrix of AMEs where element (k,j) is the average marginal effect of variable k on P(y = j).

Formula: dP(y=j)/dxk = pj * (beta{j,k} - summ pm * beta{m,k}) where beta for the base category = 0.

Key property: AMEs sum to zero across categories for each variable.

Returns

Named tuple with fields:

• effects::Matrix{T} – K x J matrix of AMEs
• varnames::Vector{String} – variable names
• categories::Vector – category labels

References

• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press, ch. 15.
• Train, K. E. (2009). Discrete Choice Methods with Simulation. 2nd ed. Cambridge Univ. Press.

marginal_effects(m::PanelLogitModel{T}; conf_level=0.95, n_quadrature=12) -> MarginalEffects{T}

Compute average marginal effects from a panel logit model with delta-method SEs.

Formulas by model type

• Pooled: AMEj = (1/n) sumi Lambda(xi'beta)(1-Lambda(xi'beta)) * beta_j
• FE (conditional): AMEj = (1/n) sumi Lambda(xi'beta)(1-Lambda(xi'beta)) * beta_j (evaluated at observed X without entity effects)
• RE: AMEj = (1/n) sumi integral Lambda(xi'beta + sigmausqrt(2)u)(1-Lambda(...)) * beta_j * w(u) du (integrated over random effect using Gauss-Hermite quadrature)
• CRE: Same as RE but only reports AMEs for original variables (not group-mean augmented ones)

Arguments

• m::PanelLogitModel{T} – estimated panel logit model
• conf_level::Real – confidence level for CIs (default 0.95)
• n_quadrature::Int – number of Gauss-Hermite quadrature points for RE/CRE (default 12)

Returns

MarginalEffects{T} with AMEs, delta-method SEs, z-stats, p-values, and CIs.

References

• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press.
• Cameron, A. C. & Trivedi, P. K. (2005). Microeconometrics. Cambridge University Press.

marginal_effects(m::PanelProbitModel{T}; conf_level=0.95, n_quadrature=12) -> MarginalEffects{T}

Compute average marginal effects from a panel probit model with delta-method SEs.

Formulas by model type

• Pooled: AMEj = (1/n) sumi phi(xi'beta) * betaj
• RE: AMEj = (1/n) sumi integral phi(xi'beta + sigmausqrt(2)u) * beta_j * w(u) du
• CRE: Same as RE but only reports AMEs for original variables

Arguments

• m::PanelProbitModel{T} – estimated panel probit model
• conf_level::Real – confidence level for CIs (default 0.95)
• n_quadrature::Int – Gauss-Hermite quadrature points for RE/CRE (default 12)

Returns

MarginalEffects{T} with AMEs, delta-method SEs, z-stats, p-values, and CIs.

odds_ratio(m::LogitModel{T}; conf_level=0.95) -> NamedTuple

Compute odds ratios from a logistic regression model with delta-method CIs.

ORj = exp(betaj). Standard error via delta method: SE(OR) = OR * SE(beta). Confidence intervals on log scale: exp(beta +/- z * SE(beta)).

Arguments

• m::LogitModel{T} — estimated logit model
• conf_level::Real — confidence level (default 0.95)

Returns

Named tuple with fields:

• or::Vector{T} — odds ratios exp(beta)
• se::Vector{T} — delta-method standard errors
• ci_lower::Vector{T} — lower CI bounds
• ci_upper::Vector{T} — upper CI bounds
• varnames::Vector{String} — variable names

Examples

m = estimate_logit(y, X)
result = odds_ratio(m)
result.or       # odds ratios
result.ci_lower # lower 95% CI

References

• Agresti, A. (2002). Categorical Data Analysis. 2nd ed. Wiley, ch. 5.

vif(m::RegModel{T}) -> Vector{T}

Compute Variance Inflation Factors for each non-intercept regressor.

VIFj = 1 / (1 - R^2j), where R^2j is the R-squared from regressing xj on all other regressors (excluding the intercept column). A VIF > 10 indicates severe multicollinearity (Belsley, Kuh & Welsch 1980).

Arguments

• m::RegModel{T} — estimated OLS/WLS model

Returns

Vector{T} of VIF values, one per non-intercept regressor. The order matches m.varnames with the intercept column removed.

Examples

m = estimate_reg(y, X; varnames=["const", "x1", "x2"])
v = vif(m)  # VIF for x1, x2

References

• Belsley, D. A., Kuh, E. & Welsch, R. E. (1980). Regression Diagnostics. Wiley.
• Greene, W. H. (2018). Econometric Analysis. 8th ed. Pearson, ch. 4.

classification_table(m::Union{LogitModel{T},ProbitModel{T}}; threshold=0.5) -> Dict{String,Any}

Compute a classification table (confusion matrix) and summary metrics for a binary response model.

Arguments

• m — estimated LogitModel or ProbitModel
• threshold::Real — classification threshold (default 0.5)

Returns

Dict with keys:

• "confusion" — 2x2 confusion matrix [[TN, FP], [FN, TP]]
• "accuracy" — (TP + TN) / N
• "sensitivity" — TP / (TP + FN) (true positive rate / recall)
• "specificity" — TN / (TN + FP) (true negative rate)
• "precision" — TP / (TP + FP) (positive predictive value)
• "f1_score" — 2 * precision * recall / (precision + recall)
• "n" — number of observations
• "threshold" — classification threshold used

Examples

m = estimate_logit(y, X)
ct = classification_table(m)
ct["accuracy"]    # overall accuracy
ct["confusion"]   # 2x2 confusion matrix

References

• Agresti, A. (2002). Categorical Data Analysis. 2nd ed. Wiley.

estimate_ologit(y, X; cov_type=:ols, varnames=nothing, clusters=nothing, maxiter=200, tol=1e-8) -> OrderedLogitModel{T}

Estimate an ordered logistic regression model via maximum likelihood (Newton-Raphson).

Model

Cumulative link model: P(y <= j | x) = Logistic(alpha_j - x' beta). X should NOT include an intercept column – it is absorbed into cutpoints.

Arguments

• y::AbstractVector – ordinal dependent variable (will be remapped to 1:J)
• X::AbstractMatrix{T} – regressor matrix (n x K, no intercept)
• cov_type::Symbol – covariance estimator: :ols (MLE), :hc1 (sandwich)
• varnames::Union{Nothing,Vector{String}} – coefficient names (auto-generated if nothing)
• clusters::Union{Nothing,AbstractVector} – cluster assignments (for :cluster)
• maxiter::Int – maximum Newton-Raphson iterations (default 200)
• tol – convergence tolerance (default 1e-8)

Returns

OrderedLogitModel{T} with estimated coefficients, cutpoints, and joint vcov.

Examples

using MacroEconometricModels, Random, Distributions
rng = MersenneTwister(42)
n = 1000
X = randn(rng, n, 2)
xb = X * [1.0, -0.5]
p = 1 ./ (1 .+ exp.(-([0.0 1.5]' .- xb)))
u = rand(rng, n)
y = [u[i] < p[1,i] ? 1 : u[i] < p[2,i] ? 2 : 3 for i in 1:n]
m = estimate_ologit(y, X; varnames=["x1", "x2"])
report(m)

References

• McCullagh, P. (1980). JRSS B 42(2), 109-142.
• Agresti, A. (2010). Analysis of Ordinal Categorical Data. 2nd ed. Wiley.

estimate_oprobit(y, X; cov_type=:ols, varnames=nothing, clusters=nothing, maxiter=200, tol=1e-8) -> OrderedProbitModel{T}

Estimate an ordered probit regression model via maximum likelihood (Newton-Raphson).

Model

Cumulative link model: P(y <= j | x) = Phi(alpha_j - x' beta). X should NOT include an intercept column – it is absorbed into cutpoints.

Arguments

Same as estimate_ologit.

Returns

OrderedProbitModel{T} with estimated coefficients, cutpoints, and joint vcov.

Examples

using MacroEconometricModels, Random, Distributions
rng = MersenneTwister(42)
n = 1000
X = randn(rng, n, 2)
xb = X * [0.8, -0.5]
d = Normal()
p = cdf.(d, [0.0 1.0]' .- xb)
u = rand(rng, n)
y = [u[i] < p[1,i] ? 1 : u[i] < p[2,i] ? 2 : 3 for i in 1:n]
m = estimate_oprobit(y, X; varnames=["x1", "x2"])
report(m)

References

• McCullagh, P. (1980). JRSS B 42(2), 109-142.
• Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data. 2nd ed. MIT Press.

estimate_mlogit(y, X; cov_type=:ols, varnames=nothing, clusters=nothing, maxiter=200, tol=1e-8) -> MultinomialLogitModel{T}

Estimate a multinomial logistic regression model via maximum likelihood (Newton-Raphson).

Model

Multinomial logit: P(y = j | x) = exp(x' betaj) / sumk exp(x' betak), with beta1 = 0 (base category).

Arguments

• y::AbstractVector – categorical dependent variable (will be remapped to 1:J)
• X::AbstractMatrix{T} – regressor matrix (n x K)
• cov_type::Symbol – covariance estimator: :ols (MLE), :hc0, :hc1 (sandwich), :cluster
• varnames::Union{Nothing,Vector{String}} – coefficient names (auto-generated if nothing)
• clusters::Union{Nothing,AbstractVector} – cluster assignments (for :cluster)
• maxiter::Int – maximum Newton-Raphson iterations (default 200)
• tol – convergence tolerance (default 1e-8)

Returns

MultinomialLogitModel{T} with estimated coefficients and joint vcov.

Examples

using MacroEconometricModels, Random
rng = MersenneTwister(42)
n = 1000
X = [ones(n) randn(rng, n, 2)]
beta_true = [0.5 -0.3; 1.0 -0.5; -0.5 0.8]  # K=3 x (J-1)=2
V = X * beta_true
P = exp.(V) ./ (1 .+ sum(exp.(V), dims=2))
# ... generate y from probabilities
m = estimate_mlogit(y, X; varnames=["const", "x1", "x2"])
report(m)

References

• McFadden, D. (1974). Conditional logit analysis of qualitative choice behavior.
• Train, K. E. (2009). Discrete Choice Methods with Simulation. 2nd ed. Cambridge Univ. Press.

brant_test(m::OrderedLogitModel{T}) -> NamedTuple

Brant test of the proportional odds (parallel regression) assumption for an ordered logit model.

For each cutpoint j = 1,...,J-1, fits a binary logit (y <= j vs y > j). Under H0 (proportional odds), all binary logit coefficients should be equal.

Algorithm

1. Fit J-1 binary logits splitting at each cutpoint
2. Compute Wald statistic comparing each binary logit's beta to the pooled ordered logit estimate
3. Overall test: chi-squared with K*(J-2) degrees of freedom
4. Per-variable tests: chi-squared with (J-2) df each

Returns

Named tuple with fields:

• statistic::T – overall Wald test statistic
• pvalue::T – overall p-value (chi-squared)
• df::Int – degrees of freedom K*(J-2)
• per_variable::Vector{T} – per-variable p-values (length K)
• binary_coefs::Matrix{T} – K x (J-1) matrix of binary logit coefficients

References

• Brant, R. (1990). Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics 46, 1171-1178.

hausman_iia(m::MultinomialLogitModel{T}; omit_category::Int) -> NamedTuple

Hausman-McFadden test of the Independence of Irrelevant Alternatives (IIA) assumption for a multinomial logit model.

Re-estimates the model excluding observations where y = omit_category, then compares coefficients for the remaining categories with the full model using a Hausman-type statistic.

Arguments

• m::MultinomialLogitModel{T} – estimated full multinomial logit model
• omit_category::Int – category to omit (1-based index into m.categories)

Returns

Named tuple with fields:

• statistic::T – Hausman test statistic
• pvalue::T – p-value (chi-squared)
• df::Int – degrees of freedom
• omitted_category – the omitted category label

References

• Hausman, J. & McFadden, D. (1984). Specification Tests for the Multinomial Logit Model. Econometrica 52(5), 1219-1240.

estimate_var(Y::AbstractMatrix{T}, p::Int; check_stability::Bool=true) -> VARModel{T}

Estimate VAR(p) via OLS: Yₜ = c + A₁Yₜ₋₁ + ... + AₚYₜ₋ₚ + uₜ.

Arguments

• Y: Data matrix (T × n)
• p: Number of lags
• check_stability: If true (default), warns if estimated VAR is non-stationary

Returns

VARModel with estimated coefficients, residuals, covariance matrix, and information criteria.

Estimate VAR from DataFrame. Use vars to select columns.

Select optimal lag order via information criterion (:aic, :bic, :hqic).

Covariance of vectorized coefficients: Σ ⊗ (X'X)⁻¹.

StatsAPI.predict(m::RegModel{T}, X_new::AbstractMatrix) -> Vector{T}

Predict fitted values for new data from a linear regression model.

Returns X_new * beta.

Arguments

• m::RegModel{T} — estimated OLS/WLS/IV model
• X_new::AbstractMatrix — new regressor matrix (n_new x k)

Returns

Vector{T} of predicted values.

Examples

m = estimate_reg(y, X)
y_hat = predict(m, X_test)

StatsAPI.predict(m::LogitModel{T}, X_new::AbstractMatrix) -> Vector{T}

Predict probabilities for new data from a logistic regression model.

Returns 1 / (1 + exp(-X_new * beta)).

Arguments

• m::LogitModel{T} — estimated logit model
• X_new::AbstractMatrix — new regressor matrix (n_new x k)

Returns

Vector{T} of predicted probabilities in (0, 1).

Examples

m = estimate_logit(y, X)
p_hat = predict(m, X_test)

StatsAPI.predict(m::ProbitModel{T}, X_new::AbstractMatrix) -> Vector{T}

Predict probabilities for new data from a probit regression model.

Returns Phi(X_new * beta), where Phi is the standard normal CDF.

Arguments

• m::ProbitModel{T} — estimated probit model
• X_new::AbstractMatrix — new regressor matrix (n_new x k)

Returns

Vector{T} of predicted probabilities in (0, 1).

Examples

m = estimate_probit(y, X)
p_hat = predict(m, X_test)

StatsAPI.predict(m::OrderedLogitModel{T}, X_new::AbstractMatrix) -> Matrix{T}

Predict category probabilities for new data from an ordered logit model.

Returns an n_new x J probability matrix where each row sums to 1.

StatsAPI.predict(m::OrderedProbitModel{T}, X_new::AbstractMatrix) -> Matrix{T}

Predict category probabilities for new data from an ordered probit model.

Returns an n_new x J probability matrix where each row sums to 1.

StatsAPI.predict(m::MultinomialLogitModel{T}, X_new::AbstractMatrix) -> Matrix{T}

Predict category probabilities for new data from a multinomial logit model.

Returns an n_new x J probability matrix where each row sums to 1.

In-sample fitted values.

Out-of-sample forecasts for steps periods.

Predicted values: F * Λ'.

Predicted values: F * Λ'.

Predicted values (common component).

predict(model::AbstractARIMAModel, h::Int) -> Vector{T}

Return h-step ahead point forecasts (without confidence intervals).

Conditional variance series $\hat{\sigma}^2_t$.

Conditional variance series $\hat{\sigma}^2_t$.

Conditional variance series $\hat{\sigma}^2_t$.

Conditional variance series $\hat{\sigma}^2_t$.

Posterior mean volatility series $\hat{\sigma}^2_t$.

R² for each equation.

R² for each variable.

R² for each variable.

R² for each variable.

Gaussian log-likelihood.

Log-likelihood of the fitted model.

Maximized log-likelihood.

Maximized log-likelihood.

Maximized log-likelihood.

Maximized log-likelihood.

Confidence intervals for RegModel coefficients (t-distribution).

Confidence intervals for LogitModel coefficients (normal distribution).

Confidence intervals for ProbitModel coefficients (normal distribution).

Confidence intervals at given level (default 95%).

Confidence intervals for PanelRegModel coefficients (t-distribution).

Confidence intervals for PanelIVModel coefficients (t-distribution).

Confidence intervals for PanelLogitModel coefficients (normal distribution).

Confidence intervals for PanelProbitModel coefficients (normal distribution).

estimate_bvar(Y, p; n_draws=1000, sampler=:direct, burnin=0, thin=1,
              prior=:normal, hyper=nothing) -> BVARPosterior

Estimate Bayesian VAR via conjugate Normal-Inverse-Wishart posterior.

Model

Y = X B + E,    E ~ MN(0, Σ, I_T)
Prior: Σ ~ IW(ν₀, S₀),  vec(B)|Σ ~ N(b₀, Σ ⊗ V₀)

Samplers

• :direct (default) — i.i.d. draws from analytical posterior. burnin and thin are ignored.
• :gibbs — Standard two-block Gibbs sampler. burnin defaults to 200 if not specified.

Arguments

• Y::AbstractMatrix: T × n data matrix
• p::Int: Number of lags

Keyword Arguments

• n_draws::Int=1000: Number of posterior draws to keep
• sampler::Symbol=:direct: Sampling algorithm (:direct or :gibbs)
• burnin::Int=0: Burnin period (only for :gibbs; defaults to 200 when sampler=:gibbs and burnin=0)
• thin::Int=1: Thinning interval (only for :gibbs)
• prior::Symbol=:normal: Prior type (:normal or :minnesota)
• hyper::Union{Nothing,MinnesotaHyperparameters}=nothing: Minnesota hyperparameters. When prior=:minnesota and hyper=nothing, tau is automatically optimized via marginal likelihood maximization (Giannone, Lenza & Primiceri 2015). Pass an explicit MinnesotaHyperparameters(...) to use fixed values instead.

Returns

BVARPosterior{T} containing coefficient and covariance draws.

Example

Y = randn(200, 3)
post = estimate_bvar(Y, 2; n_draws=1000)
post_mn = estimate_bvar(Y, 2; prior=:minnesota, n_draws=500)

VARModel with posterior mean parameters.

VARModel with posterior median parameters.

Compute level residual covariance Sigma for structural analysis.

Draw from Inverse-Wishart(ν, S) distribution. Uses the Bartlett decomposition: if X ~ W(ν, S⁻¹), then X⁻¹ ~ IW(ν, S).

_estimate_favar_bayesian(X, Y_key, r, p, n_draws, burnin, pvn, Y_key_indices, aug_varnames) -> BayesianFAVAR

Bayesian FAVAR via Gibbs sampling (Bernanke, Boivin & Eliasz 2005, Section IV).

Algorithm:

1. Initialize factors via PCA
2. Gibbs sampler iterates:
   • Draw (B, Σ) | F, Y_key from the VAR posterior
   • Draw Λ | F, X equation-by-equation
   • Draw F | Λ, B, Σ, X, Y_key via regression posterior
3. Store post-burnin draws

_find_key_indices(X, Y_key) -> Vector{Int}

Try to identify which columns of X correspond to Y_key by exact column matching. Returns empty vector if no match is found.

_remove_double_counting(F_raw, Y_key) -> Matrix{T}

Remove the component of Fraw that is spanned by Ykey (BBE 2005 Step 2).

For each factor column, regress it on Y_key via OLS and keep the residual. This yields "slow-moving" factors that are orthogonal to the key variables.

Draw i.i.d. samples from the conjugate Normal-Inverse-Wishart posterior.

Two-block Gibbs sampler for Normal-Inverse-Wishart posterior.

Windmeijer (2005) finite-sample correction for equation eq.

estimate_bvar(Y, p; n_draws=1000, sampler=:direct, burnin=0, thin=1,
              prior=:normal, hyper=nothing) -> BVARPosterior

Estimate Bayesian VAR via conjugate Normal-Inverse-Wishart posterior.

Model

Y = X B + E,    E ~ MN(0, Σ, I_T)
Prior: Σ ~ IW(ν₀, S₀),  vec(B)|Σ ~ N(b₀, Σ ⊗ V₀)

Samplers

• :direct (default) — i.i.d. draws from analytical posterior. burnin and thin are ignored.
• :gibbs — Standard two-block Gibbs sampler. burnin defaults to 200 if not specified.

Arguments

• Y::AbstractMatrix: T × n data matrix
• p::Int: Number of lags

Keyword Arguments

• n_draws::Int=1000: Number of posterior draws to keep
• sampler::Symbol=:direct: Sampling algorithm (:direct or :gibbs)
• burnin::Int=0: Burnin period (only for :gibbs; defaults to 200 when sampler=:gibbs and burnin=0)
• thin::Int=1: Thinning interval (only for :gibbs)
• prior::Symbol=:normal: Prior type (:normal or :minnesota)
• hyper::Union{Nothing,MinnesotaHyperparameters}=nothing: Minnesota hyperparameters. When prior=:minnesota and hyper=nothing, tau is automatically optimized via marginal likelihood maximization (Giannone, Lenza & Primiceri 2015). Pass an explicit MinnesotaHyperparameters(...) to use fixed values instead.

Returns

BVARPosterior{T} containing coefficient and covariance draws.

Example

Y = randn(200, 3)
post = estimate_bvar(Y, 2; n_draws=1000)
post_mn = estimate_bvar(Y, 2; prior=:minnesota, n_draws=500)

estimate_favar(X, Y_key, r, p; method=:two_step, panel_varnames=nothing) -> FAVARModel
estimate_favar(X, key_indices::Vector{Int}, r, p; ...) -> FAVARModel

Estimate a Factor-Augmented VAR (Bernanke, Boivin & Eliasz, 2005).

Two-Step Algorithm (BBE 2005)

1. Extract r factors from panel X via PCA (estimate_factors)
2. Remove double-counting: regress each factor on Y_key, use residuals as slow-moving factors F_tilde
3. Estimate VAR(p) on augmented system [F_tilde, Y_key]

Arguments

• X: Panel data matrix (T x N) — large cross-section of macroeconomic variables
• Y_key: Key observed variables (T x n_key) or column indices into X
• r: Number of latent factors to extract
• p: VAR lag order

Keyword Arguments

• method::Symbol=:two_step: Estimation method (currently only :two_step)
• panel_varnames::Union{Nothing, Vector{String}}=nothing: Names for the N panel variables
• n_draws::Int=5000: Number of posterior draws (for Bayesian, reserved for future use)
• burnin::Int=1000: Burn-in draws (for Bayesian, reserved for future use)

Returns

FAVARModel{T} with estimated VAR on [Ftilde, Ykey], factor loadings, and panel variable metadata for panel-wide IRF mapping.

Example

X = randn(200, 50)    # 200 obs, 50 variables
Y_key = X[:, [1, 5]]  # 2 key variables
favar = estimate_favar(X, Y_key, 3, 2)

# Or using column indices:
favar = estimate_favar(X, [1, 5], 3, 2)

# IRFs (via to_var delegation):
irf_result = irf(favar, 20)

# Panel-wide IRFs (all 50 variables):
panel_irfs = favar_panel_irf(favar, irf_result)

estimate_favar(X, key_indices::Vector{Int}, r, p; kwargs...) -> FAVARModel or BayesianFAVAR

Estimate FAVAR where key variables are specified as column indices of X.

estimate_pvar(d::PanelData{T}, p::Int; kwargs...) -> PVARModel{T}

Estimate Panel VAR(p) via GMM.

Arguments

• d::PanelData{T} — balanced panel data (use xtset to construct)
• p::Int — number of lags

Keyword Arguments

• dependent_vars::Union{Vector{String},Nothing}=nothing — endogenous variable names (default: all)
• predet_vars::Vector{String}=String[] — predetermined variable names
• exog_vars::Vector{String}=String[] — strictly exogenous variable names
• transformation::Symbol=:fd — :fd (first-difference) or :fod (forward orthogonal deviations)
• steps::Symbol=:twostep — :onestep, :twostep, or :mstep (iterated)
• system_instruments::Bool=false — if true, use System GMM (Blundell-Bond)
• system_constant::Bool=true — include constant in level equation (System GMM)
• min_lag_endo::Int=2 — minimum instrument lag for endogenous variables
• max_lag_endo::Int=99 — maximum instrument lag (99 = all available)
• collapse::Bool=false — collapse instruments to limit proliferation
• pca_instruments::Bool=false — apply PCA reduction to instruments
• pca_max_components::Int=0 — max PCA components (0 = auto)
• max_iter::Int=100 — max iterations for iterated GMM

Returns

• PVARModel{T} with coefficient estimates, robust standard errors, and GMM internals

References

• Arellano, M. & Bond, S. (1991). Review of Economic Studies 58(2), 277-297.
• Blundell, R. & Bond, S. (1998). Journal of Econometrics 87(1), 115-143.
• Windmeijer, F. (2005). Journal of Econometrics 126(1), 25-51.

Examples

pd = xtset(df, :id, :time)
model = estimate_pvar(pd, 2; steps=:twostep)
model = estimate_pvar(pd, 1; system_instruments=true, steps=:twostep)

estimate_pvar_feols(d::PanelData{T}, p::Int; kwargs...) -> PVARModel{T}

Estimate Panel VAR(p) via Fixed Effects OLS (within estimator).

Simpler than GMM but subject to Nickell (1981) bias when T is small relative to N. Uses within-group demeaning to remove fixed effects, then pooled OLS with cluster-robust standard errors.

Arguments

• d::PanelData{T} — panel data
• p::Int — number of lags

Keyword Arguments

• dependent_vars::Union{Vector{String},Nothing}=nothing — endogenous variable names
• predet_vars::Vector{String}=String[] — predetermined variable names
• exog_vars::Vector{String}=String[] — strictly exogenous variable names

Examples

pd = xtset(df, :id, :time)
model = estimate_pvar_feols(pd, 2)

Estimate VAR from DataFrame. Use vars to select columns.

estimate_var(Y::AbstractMatrix{T}, p::Int; check_stability::Bool=true) -> VARModel{T}

Estimate VAR(p) via OLS: Yₜ = c + A₁Yₜ₋₁ + ... + AₚYₜ₋ₚ + uₜ.

Arguments

• Y: Data matrix (T × n)
• p: Number of lags
• check_stability: If true (default), warns if estimated VAR is non-stationary

Returns

VARModel with estimated coefficients, residuals, covariance matrix, and information criteria.

extract_chain_parameters(post::BVARPosterior) -> (b_vecs, sigmas)

Extract coefficient vectors and covariance matrices from posterior draws. Returns matrices compatible with downstream parameters_to_model calls.

• b_vecs: n_draws × (k*n) matrix of vectorized coefficients
• sigmas: n_draws × (n*n) matrix of vectorized covariance matrices

forecast(post::BVARPosterior, h; reps=nothing, conf_level=0.95, point_estimate=:mean) -> BVARForecast{T}

Forecast from Bayesian VAR using posterior draws.

For each posterior draw, iterates the VAR recursion forward h steps. Returns posterior mean forecast with credible intervals from the distribution of forecast paths across draws.

Arguments

• post: BVAR posterior
• h: Forecast horizon

Keyword Arguments

• reps: Number of posterior draws to use (default: all)
• conf_level: Credible interval level (default: 0.95)
• point_estimate: :mean (default) or :median for central tendency

Returns

BVARForecast{T} with posterior mean forecast and credible intervals.

Example

post = estimate_bvar(Y, 4; n_draws=1000)
fc = forecast(post, 12)

forecast(model::VARModel, h; ci_method=:bootstrap, reps=500, conf_level=0.95) -> VARForecast{T}

Forecast from VAR model for h steps ahead with optional bootstrap CIs.

Iterates the VAR recursion forward using the estimated coefficients and the last p observations as initial conditions.

Arguments

• model: Estimated VAR model
• h: Forecast horizon (number of steps ahead)
• ci_method: :bootstrap (default) or :none
• reps: Number of bootstrap replications (default 500)
• conf_level: Confidence level for intervals (default 0.95)

Returns

VARForecast{T} with point forecasts and CIs.

Example

model = estimate_var(Y, 4)
fc = forecast(model, 12)  # 12-step ahead forecast with bootstrap CIs

Convert chain parameters to VARModel. Provide data for residual computation.

VARModel with posterior mean parameters.

VARModel with posterior median parameters.

Select optimal lag order via information criterion (:aic, :bic, :hqic).

Confidence intervals at given level (default 95%).

Gaussian log-likelihood.

In-sample fitted values.

Out-of-sample forecasts for steps periods.

R² for each equation.

Covariance of vectorized coefficients: Σ ⊗ (X'X)⁻¹.

point_forecast(f::AbstractForecastResult)

Return the point forecast values (Vector or Matrix).

Most subtypes store the point forecast in a .forecast field; this default accessor returns that field. Overrides exist for VECMForecast (.levels) and FactorForecast (.observables).

forecast_horizon(f::AbstractForecastResult)

Return the forecast horizon (number of steps ahead).

lower_bound(f::AbstractForecastResult)

Return the lower confidence interval bound.

Most subtypes store this in .ci_lower; an override exists for FactorForecast (.observables_lower).

upper_bound(f::AbstractForecastResult)

Return the upper confidence interval bound.

Most subtypes store this in .ci_upper; an override exists for FactorForecast (.observables_upper).

gen_dummy_obs(Y, p, hyper) -> (Y_dummy, X_dummy)

Generate Minnesota prior dummy observations. Hyperparameters: tau (tightness), decay, lambda (sum-of-coef), mu (co-persistence), omega.

log_marginal_likelihood(Y, p, hyper) -> T

Closed-form log marginal likelihood for BVAR with Minnesota prior.

Optimize tau via grid search on marginal likelihood.

Full grid search over tau, lambda, mu.

estimate_vecm(Y, p; rank=:auto, deterministic=:constant, method=:johansen, significance=0.05)

Estimate a Vector Error Correction Model.

Arguments

• Y: Data matrix (T × n) in levels
• p: Underlying VAR order (VECM has p-1 lagged differences)
• rank: Cointegrating rank; :auto (default) selects via Johansen trace test, or specify an integer
• deterministic: :none, :constant (default), or :trend
• method: :johansen (default) or :engle_granger (bivariate, rank=1 only)
• significance: Significance level for automatic rank selection (default 0.05)

Returns

VECMModel with estimated α, β, Γ matrices, residuals, and diagnostics.

Example

Y = cumsum(randn(200, 3), dims=1)
Y[:, 2] = Y[:, 1] + 0.1 * randn(200)
m = estimate_vecm(Y, 2)

to_var(vecm::VECMModel) -> VARModel

Convert VECM to VAR in levels representation.

The VECM: ΔYₜ = ΠYₜ₋₁ + Σᵢ ΓᵢΔYₜ₋ᵢ + μ + uₜ maps to VAR: Yₜ = c + A₁Yₜ₋₁ + ... + AₚYₜ₋ₚ + uₜ

with:

• A₁ = Π + Iₙ + Γ₁
• Aᵢ = Γᵢ - Γᵢ₋₁ for i = 2, ..., p-1
• Aₚ = -Γₚ₋₁

to_var(favar::FAVARModel) -> VARModel

Convert a FAVAR to a standard VARModel for use with irf, fevd, historical_decomposition, and other structural analysis methods.

This follows the same delegation pattern as to_var(vecm::VECMModel).

select_vecm_rank(Y, p; criterion=:trace, significance=0.05) -> Int

Select cointegrating rank using the Johansen trace or max-eigenvalue test.

cointegrating_rank(m::VECMModel) -> Int

Return the cointegrating rank r of the VECM.

granger_causality_vecm(vecm, cause, effect) -> VECMGrangerResult

Test Granger causality from variable cause to variable effect in a VECM.

Three tests are computed:

1. Short-run: Wald test on Γ coefficients of the cause variable in the effect equation
2. Long-run: Wald test on α coefficients in the effect equation (error correction channel)
3. Strong: Joint test of both short-run and long-run

Arguments

• vecm: Estimated VECMModel
• cause: Index of the causing variable
• effect: Index of the effect variable

Returns

VECMGrangerResult with test statistics, p-values, and degrees of freedom.

fevd(vecm::VECMModel, horizon; kwargs...) -> FEVD

Compute FEVD for a VECM by converting to VAR representation.

historical_decomposition(vecm::VECMModel, horizon; kwargs...) -> HistoricalDecomposition

Compute historical decomposition for a VECM by converting to VAR representation.

irf(vecm::VECMModel, horizon; kwargs...) -> ImpulseResponse

Compute IRFs for a VECM by converting to VAR representation. All identification methods (Cholesky, sign, narrative, etc.) are supported.

forecast(vecm::VECMModel, h; ci_method=:none, reps=500, conf_level=0.95) -> VECMForecast

Forecast from a VECM by iterating the VECM equations in levels.

Unlike VAR forecasting, this preserves the cointegrating relationships in the forecast path.

Arguments

• vecm: Estimated VECM
• h: Forecast horizon
• ci_method: :none (default), :bootstrap, or :simulation
• reps: Number of bootstrap/simulation replications (default 500)
• conf_level: Confidence level (default 0.95)

Returns

VECMForecast with level and difference forecasts, plus CIs if requested.

compute_Q(model, method, horizon, check_func, narrative_check;
          max_draws=100, transition_var=nothing, regime_indicator=nothing)

Compute identification matrix Q for structural VAR analysis.

Methods

• :cholesky — Cholesky decomposition (recursive ordering)
• :sign — Sign restrictions (requires check_func)
• :narrative — Narrative restrictions (requires check_func and narrative_check)
• :long_run — Long-run restrictions (Blanchard-Quah)
• :fastica — FastICA (Hyvärinen 1999)
• :jade — JADE (Cardoso 1999)
• :sobi — SOBI (Belouchrani et al. 1997)
• :dcov — Distance covariance ICA (Matteson & Tsay 2017)
• :hsic — HSIC independence ICA (Gretton et al. 2005)
• :student_t — Student-t ML (Lanne et al. 2017)
• :mixture_normal — Mixture of normals ML (Lanne et al. 2017)
• :pml — Pseudo-ML (Gouriéroux et al. 2017)
• :skew_normal — Skew-normal ML (Lanne & Luoto 2020)
• :nongaussian_ml — Unified non-Gaussian ML dispatcher (default: Student-t)
• :markov_switching — Markov-switching heteroskedasticity (Lütkepohl & Netšunajev 2017)
• :garch — GARCH-based heteroskedasticity (Normandin & Phaneuf 2004)
• :smooth_transition — Smooth-transition heteroskedasticity (requires transition_var)
• :external_volatility — External volatility regimes (requires regime_indicator)

Keyword Arguments

• max_draws::Int=100: Maximum draws for sign/narrative identification
• transition_var::Union{Nothing,AbstractVector}=nothing: Transition variable for :smooth_transition
• regime_indicator::Union{Nothing,AbstractVector{Int}}=nothing: Regime indicator for :external_volatility

compute_irf(model, Q, horizon) -> Array{T,3}

Compute IRFs for rotation matrix Q. Returns (horizon × n × n) array. IRF[h, i, j] = response of variable i to shock j at horizon h-1.

Compute structural shocks: εₜ = Q'L⁻¹uₜ.

Generate random orthogonal matrix via QR decomposition (Haar measure).

Identify via Cholesky decomposition (recursive ordering). Returns L where Σ = LL'.

Identify via long-run restrictions: long-run cumulative impact matrix is lower triangular.

identify_narrative(model, horizon, sign_check, narrative_check; max_draws=1000)

Combine sign and narrative restrictions. Returns (Q, irf, shocks).

identify_sign(model, horizon, check_func; max_draws=1000, store_all=false)

Find Q satisfying sign restrictions via random draws.

With store_all=false (default), returns (Q, irf) — the first valid rotation. With store_all=true, returns a SignIdentifiedSet containing ALL accepted rotations and their IRFs (Baumeister & Hamilton, 2015).

irf_bounds(s::SignIdentifiedSet{T}; quantiles=[0.16, 0.84]) -> (lower, upper)

Compute pointwise bounds (or quantile bands) over the identified set.

irf_median(s::SignIdentifiedSet{T}) -> Array{T,3}

Compute pointwise median IRF over the identified set.

identify_arias(model, restrictions, horizon; n_draws=1000, n_rotations=1000, compute_weights=true) -> AriasSVARResult

Identify SVAR using Arias et al. (2018) with zero and sign restrictions.

Uses importance sampling with draw-dependent weights (Proposition 4) for zero+sign restriction combinations. For pure sign restrictions, draws uniformly from O(n) with unit weights.

Keywords

• n_draws::Int=1000: Target number of accepted draws
• n_rotations::Int=1000: Maximum attempts per target draw
• compute_weights::Bool=true: Compute importance weights (set false for faster exploratory analysis)

identify_arias_bayesian(post::BVARPosterior, restrictions, horizon; data=nothing, n_rotations=100, quantiles=[0.16,0.5,0.84], compute_weights=true)

Apply Arias identification to each posterior draw. Returns IRF quantiles, mean, acceptance rates.

Creates the _AriasSVARSetup once (W matrices fixed across all posterior draws) for consistency.

Create zero restriction: variable doesn't respond to shock at horizon.

Create sign restriction: variable response has given sign (:positive/:negative) at horizon.

irf_bounds(s::SignIdentifiedSet{T}; quantiles=[0.16, 0.84]) -> (lower, upper)

Compute pointwise bounds (or quantile bands) over the identified set.

irf_median(s::SignIdentifiedSet{T}) -> Array{T,3}

Compute pointwise median IRF over the identified set.

Compute weighted mean IRF from AriasSVARResult.

Compute weighted IRF percentiles from AriasSVARResult.

identify_uhlig(model::VARModel{T}, restrictions::SVARRestrictions, horizon::Int;
    n_starts=50, n_refine=10, max_iter_coarse=500, max_iter_fine=2000,
    tol_coarse=1e-4, tol_fine=1e-8) -> UhligSVARResult{T}

Identify SVAR using Mountford & Uhlig (2009) penalty function approach.

Uses Nelder-Mead optimization over spherical coordinates to find the rotation matrix $Q$ that best satisfies sign restrictions, with zero restrictions enforced as hard constraints via null-space projection.

Algorithm

1. Precompute MA coefficients and Cholesky factor $L$
2. Phase 1 (coarse): n_starts Nelder-Mead runs from random $\theta_0 \in [0, 2\pi]$
3. Phase 2 (refinement): n_refine local re-optimizations from best solution
4. Build final $Q$, compute IRFs, check convergence

Keywords

• n_starts::Int=50: Number of random starting points (Phase 1)
• n_refine::Int=10: Number of local refinements (Phase 2)
• max_iter_coarse::Int=500: Max iterations per Phase 1 run
• max_iter_fine::Int=2000: Max iterations per Phase 2 run
• tol_coarse::T=1e-4: Convergence tolerance for Phase 1
• tol_fine::T=1e-8: Convergence tolerance for Phase 2

Returns

UhligSVARResult{T} with optimal rotation matrix, IRFs, penalty values, and convergence indicator.

Example

model = estimate_var(Y, 2)
restrictions = SVARRestrictions(3;
    zeros = [zero_restriction(3, 1)],
    signs = [sign_restriction(1, 1, :positive),
             sign_restriction(2, 1, :positive)]
)
result = identify_uhlig(model, restrictions, 20)

Reference: Mountford & Uhlig (2009)

Simulate IRFs for confidence intervals (bootstrap or asymptotic).

Simulate VAR data from initial conditions and innovations.

cumulative_irf(irf_result::BayesianImpulseResponse{T}) -> BayesianImpulseResponse{T}

Compute cumulative Bayesian impulse response: Σₛ₌₀ʰ IRF_s.

When raw posterior draws are available, cumulates each draw first then extracts quantiles — the statistically correct approach.

cumulative_irf(irf_result::ImpulseResponse{T}) -> ImpulseResponse{T}

Compute cumulative impulse response for VAR models: Σₛ₌₀ʰ IRF_s.

When raw bootstrap/simulation draws are available, cumulates each draw first then extracts quantiles — the statistically correct approach since quantiles are NOT additive: Qα(A+B) ≠ Qα(A) + Q_α(B).

cumulative_irf(irf::LPImpulseResponse{T}) -> LPImpulseResponse{T}

Compute cumulative impulse response: Σₛ₌₀ʰ β_s.

irf(post::BVARPosterior, horizon; method=:cholesky, quantiles=[0.16, 0.5, 0.84], point_estimate=:mean, ...)

Compute Bayesian IRFs from posterior draws with posterior quantiles. Uses posterior mean as central tendency by default (pass point_estimate=:median for median).

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

Uses process_posterior_samples and compute_posterior_quantiles from bayesian_utils.jl.

irf(slp::StructuralLP) -> ImpulseResponse

Extract the impulse response object from a structural LP result.

irf(model, horizon; method=:cholesky, ci_type=:none, reps=200, conf_level=0.95, ...)

Compute IRFs with optional confidence intervals.

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

CI types

:none, :bootstrap, :theoretical

lp_irf(Y::AbstractMatrix, shock_var::Int, horizon::Int; kwargs...) -> LPImpulseResponse

Convenience function: estimate LP and extract IRF in one call.

lp_irf(model::LPModel{T}; conf_level::Real=0.95) -> LPImpulseResponse{T}

Extract impulse response function with confidence intervals from LP model.

Compute FEVD from IRF array: decomposition[i,j,h] = cumulative MSE contribution.

fevd(model, horizon; method=:cholesky, ...) -> FEVD

Compute FEVD showing proportion of h-step forecast error variance attributable to each shock.

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

Compute historical decomposition contributions from structural shocks and MA coefficients. HD[t, i, j] = Σ{s=0}^{t-1} Θs[i, j] * ε_j(t-s)

Compute initial conditions as residual: actual - total shock contributions.

Compute structural MA coefficients Θs = Φs * P for s = 0, ..., horizon-1. Returns Vector{Matrix{T}} of length horizon.

Uses _compute_ma_coefficients from identification.jl to avoid code duplication.

contribution(hd::BayesianHistoricalDecomposition, var, shock; stat=:mean) -> Vector

Get contribution time series for specific variable and shock (Bayesian).

Arguments

• hd: Bayesian historical decomposition result
• var: Variable index (Int) or name (String)
• shock: Shock index (Int) or name (String)
• stat: :mean or quantile index (Int)

contribution(hd::HistoricalDecomposition, var, shock) -> Vector

Get contribution time series for specific variable and shock.

Arguments

• hd: Historical decomposition result
• var: Variable index (Int) or name (String)
• shock: Shock index (Int) or name (String)

Example

contrib_y1_s1 = contribution(hd, 1, 1)  # Contribution of shock 1 to variable 1
contrib_y1_s1 = contribution(hd, "Var 1", "Shock 1")

historical_decomposition(post::BVARPosterior, horizon; data=..., ...) -> BayesianHistoricalDecomposition

Compute Bayesian historical decomposition from posterior draws with posterior quantiles.

Arguments

• post::BVARPosterior: Posterior draws from estimate_bvar
• horizon::Int: Maximum horizon for MA coefficients

Keyword Arguments

• data::AbstractMatrix: Override data matrix (defaults to post.data)
• method::Symbol=:cholesky: Identification method
• quantiles::Vector{<:Real}=[0.16, 0.5, 0.84]: Posterior quantile levels
• check_func=nothing: Sign restriction check function
• narrative_check=nothing: Narrative restriction check function
• transition_var=nothing: Transition variable (for method=:smooth_transition)
• regime_indicator=nothing: Regime indicator (for method=:external_volatility)

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

Returns

BayesianHistoricalDecomposition with posterior quantiles and means.

Example

post = estimate_bvar(Y, 2; n_draws=500)
hd = historical_decomposition(post, 198)

historical_decomposition(slp::StructuralLP{T}, T_hd::Int) -> HistoricalDecomposition{T}

Compute historical decomposition from structural LP.

Uses LP-estimated IRFs as the structural MA coefficients Θ_h and the structural shocks from the underlying VAR identification.

Arguments

• slp: Structural LP result
• T_hd: Number of time periods for decomposition (≤ T_eff of underlying VAR)

Returns

HistoricalDecomposition{T} with contributions, initial conditions, and actual data.

historical_decomposition(model::VARModel, restrictions::SVARRestrictions, horizon; ...) -> BayesianHistoricalDecomposition

Compute historical decomposition using Arias et al. (2018) identification with importance weights.

Arguments

• model::VARModel: Estimated VAR model
• restrictions::SVARRestrictions: Zero and sign restrictions
• horizon::Int: Maximum horizon for MA coefficients

Keyword Arguments

• n_draws::Int=1000: Number of accepted draws
• n_rotations::Int=1000: Maximum rotation attempts per draw
• quantiles::Vector{<:Real}=[0.16, 0.5, 0.84]: Quantile levels for weighted quantiles

Returns

BayesianHistoricalDecomposition with weighted posterior quantiles and means.

Example

r = SVARRestrictions(3; signs=[sign_restriction(1, 1, :positive)])
hd = historical_decomposition(model, r, 198; n_draws=500)

historical_decomposition(model::VARModel, horizon; method=:cholesky, ...) -> HistoricalDecomposition

Compute historical decomposition for a VAR model.

Decomposes observed data into contributions from each structural shock plus initial conditions.

Arguments

• model::VARModel: Estimated VAR model
• horizon::Int: Maximum horizon for MA coefficient computation (typically T_eff)

Keyword Arguments

• method::Symbol=:cholesky: Identification method
• check_func=nothing: Sign restriction check function (for method=:sign or :narrative)
• narrative_check=nothing: Narrative restriction check function (for method=:narrative)
• max_draws::Int=1000: Maximum draws for sign/narrative identification
• transition_var=nothing: Transition variable (for method=:smooth_transition)
• regime_indicator=nothing: Regime indicator (for method=:external_volatility)

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

Returns

HistoricalDecomposition containing:

• contributions: Shock contributions (Teff × nvars × n_shocks)
• initial_conditions: Initial condition effects (Teff × nvars)
• actual: Actual data values
• shocks: Structural shocks

Example

model = estimate_var(Y, 2)
hd = historical_decomposition(model, size(Y, 1) - 2)
verify_decomposition(hd)  # Check decomposition identity

total_shock_contribution(hd::AbstractHistoricalDecomposition, var) -> Vector

Get total contribution from all shocks to a variable over time.

verify_decomposition(hd::BayesianHistoricalDecomposition; tol=1e-6) -> Bool

Verify that mean contributions + mean initial_conditions ≈ actual (approximately, due to averaging).

verify_decomposition(hd::HistoricalDecomposition; tol=1e-10) -> Bool

Verify that contributions + initial_conditions ≈ actual.

Example

hd = historical_decomposition(model, horizon)
@assert verify_decomposition(hd) "Decomposition identity failed"

has_uncertainty(result::AbstractAnalysisResult) -> Bool

Check if the result includes uncertainty quantification (confidence intervals or posterior quantiles).

point_estimate(result::AbstractAnalysisResult)

Get the point estimate from an analysis result.

Returns the main values/estimates (IRF values, FEVD proportions, HD contributions).

report(post::BVARPosterior)

Print comprehensive Bayesian VAR posterior summary.

report(r::BaxterKingResult)

Print Baxter-King band-pass filter summary.

report(r::BeveridgeNelsonResult)

Print Beveridge-Nelson decomposition summary.

report(r::BoostedHPResult)

Print boosted HP filter summary.

report(f::FEVD)
report(f::BayesianFEVD)

Print FEVD summary with decomposition at selected horizons.

report(r::HPFilterResult)

Print HP filter summary.

report(r::HamiltonFilterResult)

Print Hamilton filter summary.

report(hd::HistoricalDecomposition)
report(hd::BayesianHistoricalDecomposition)

Print HD summary with contribution statistics.

report(irf::ImpulseResponse)
report(irf::BayesianImpulseResponse)

Print IRF summary with values at selected horizons.

report(f::VECMForecast)

Print VECM forecast summary.

report(g::VECMGrangerResult)

Print VECM Granger causality test results.

report(model::VARModel)

Print comprehensive VAR model summary including specification, per-equation coefficient estimates with standard errors and significance, information criteria, residual covariance, and stationarity check.

report(vecm::VECMModel)

Print comprehensive VECM summary including cointegrating vectors, adjustment coefficients, short-run dynamics, and diagnostics.

uncertainty_bounds(result::AbstractAnalysisResult) -> Union{Nothing, Tuple}

Get uncertainty bounds (lower, upper) if available, otherwise nothing.

Generate a block bootstrap sample from matrix Y.

Extract 3D IRF, SE, and analytical CI arrays from per-shock LP models.

_lp_predict_at_horizon(Y, response_vars, residuals_h, T_eff_h, h)

Reconstruct fitted values at horizon h via Yh − residualsh identity.

_lp_robust_vcov(X_h, U_h, cov_estimator, h)

Horizon-aware robust variance-covariance for LP regressions.

LP residuals at horizon h have MA(h-1) serial correlation (Jordà 2005), so the Newey-West bandwidth must be at least h+1. When automatic bandwidth selection (bandwidth == 0) yields a smaller value, this function enforces the floor max(m̂_NW, h+1).

For non-NW estimators (White, Driscoll-Kraay), falls through to the standard compute_block_robust_vcov.

Block bootstrap for structural LP confidence intervals.

build_control_columns!(X_h::AbstractMatrix{T}, Y::AbstractMatrix{T},
                       t_start::Int, t_end::Int, lags::Int, start_col::Int) where T

Fill control (lagged Y) columns into regressor matrix X_h. Returns the next available column index.

build_response_matrix(Y::AbstractMatrix{T}, h::Int, t_start::Int, t_end::Int,
                      response_vars::Vector{Int}) where T

Build response matrix Y_h at horizon h.

compare_var_lp(Y::AbstractMatrix{T}, horizon::Int; lags::Int=4) where T

Compare VAR-based and LP-based impulse responses.

compute_block_robust_vcov(X::AbstractMatrix{T}, U::AbstractMatrix{T},
                          cov_estimator::AbstractCovarianceEstimator) where T

Compute block-diagonal robust covariance for multi-equation system.

compute_horizon_bounds(T_obs::Int, h::Int, lags::Int) -> (t_start, t_end)

Compute valid observation bounds for horizon h.

construct_lp_matrices(Y::AbstractMatrix{T}, shock_var::Int, h::Int, lags::Int;
                      response_vars::Vector{Int}=collect(1:size(Y,2))) where T

Construct regressor and response matrices for LP regression at horizon h.

Returns: (Yh, Xh, valid_idx)

create_cov_estimator(cov_type::Symbol, ::Type{T}; bandwidth::Int=0) where T

Create covariance estimator from symbol specification using the global registry.

estimate_lp(Y::AbstractMatrix{T}, shock_var::Int, horizon::Int;
            lags::Int=4, response_vars::Vector{Int}=collect(1:size(Y,2)),
            cov_type::Symbol=:newey_west, bandwidth::Int=0,
            conf_level::Real=0.95) -> LPModel{T}

Estimate Local Projection impulse response functions (Jordà 2005).

The LP regression for horizon h: y{t+h} = αh + βh * shockt + Γh * controlst + ε_{t+h}

estimate_lp_cholesky(Y::AbstractMatrix{T}, horizon::Int;
                     lags::Int=4, cov_type::Symbol=:newey_west, kwargs...) -> Vector{LPModel{T}}

Estimate LP with Cholesky-orthogonalized shocks.

estimate_lp_multi(Y::AbstractMatrix{T}, shock_vars::Vector{Int}, horizon::Int;
                  kwargs...) -> Vector{LPModel{T}}

Estimate LP for multiple shock variables.

extract_shock_irf(B::Vector{Matrix{T}}, vcov::Vector{Matrix{T}},
                  response_vars::Vector{Int}, shock_coef_idx::Int;
                  conf_level::Real=0.95) where T

Generic IRF extraction from coefficient and covariance vectors. Works for LPModel, LPIVModel, PropensityLPModel.

structural_lp(Y::AbstractMatrix{T}, horizon::Int;
              method=:cholesky, lags=4, var_lags=nothing,
              cov_type=:newey_west, conf_level=0.95,
              ci_type=:none, reps=200,
              check_func=nothing, narrative_check=nothing,
              max_draws=1000,
              transition_var=nothing, regime_indicator=nothing) -> StructuralLP{T}

Estimate structural LP impulse responses using VAR-based identification with LP estimation.

Algorithm:

1. Estimate VAR(p) → obtain Σ and residuals
2. Compute rotation matrix Q via chosen identification method
3. Compute structural shocks εt = Q'L⁻¹ut
4. For each shock j, run LP with εj as regressor and Yeff as responses
5. Stack into 3D IRF array: irfs[h, i, j]

Arguments

• Y: T × n data matrix
• horizon: Maximum IRF horizon H

Keyword Arguments

• method: Identification method
• lags: Number of LP control lags (default: 4)
• var_lags: VAR lag order for identification (default: same as lags)
• cov_type: HAC estimator (:newey_west, :white)
• conf_level: Confidence level for CIs (default: 0.95)
• ci_type: CI method (:none, :bootstrap)
• reps: Number of bootstrap replications
• check_func: Sign restriction check function (for :sign/:narrative)
• narrative_check: Narrative check function (for :narrative)
• max_draws: Maximum draws for sign/narrative identification
• transition_var: Transition variable (for :smooth_transition)
• regime_indicator: Regime indicator (for :external_volatility)

Methods

:cholesky, :sign, :narrative, :long_run, :fastica, :jade, :sobi, :dcov, :hsic, :student_t, :mixture_normal, :pml, :skew_normal, :nongaussian_ml, :markov_switching, :garch, :smooth_transition, :external_volatility

Note: :smooth_transition requires transition_var kwarg. :external_volatility requires regime_indicator kwarg.

Returns

StructuralLP{T} with 3D IRFs, structural shocks, VAR model, and individual LP models.

References

Plagborg-Møller, M. & Wolf, C. K. (2021). "Local Projections and VARs Estimate the Same Impulse Responses." Econometrica, 89(2), 955–980.

estimate_lp_iv(Y::AbstractMatrix{T}, shock_var::Int, instruments::AbstractMatrix{T},
               horizon::Int; lags::Int=4, response_vars::Vector{Int}=collect(1:size(Y,2)),
               cov_type::Symbol=:newey_west, bandwidth::Int=0) -> LPIVModel{T}

Estimate LP with Instrumental Variables (Stock & Watson 2018) using 2SLS.

first_stage_regression(endog::AbstractVector{T}, instruments::AbstractMatrix{T},
                       controls::AbstractMatrix{T}; h::Int=0) -> NamedTuple

First-stage regression for 2SLS with HAC-robust F-statistic for instrument relevance.

At LP horizon h > 0, residuals have MA(h-1) autocorrelation by construction (Jordà 2005), so the F-statistic uses Newey-West HAC variance with bandwidth ≥ h+1.

lp_iv_irf(model::LPIVModel{T}; conf_level::Real=0.95) -> LPImpulseResponse{T}

Extract IRF from LP-IV model.

sargan_test(model::LPIVModel{T}, h::Int) -> NamedTuple

Sargan-Hansen J-test for overidentification at horizon h.

tsls_regression(Y::AbstractMatrix{T}, endog::AbstractVector{T},
                endog_fitted::AbstractVector{T}, controls::AbstractMatrix{T};
                cov_estimator::AbstractCovarianceEstimator=NeweyWestEstimator()) -> NamedTuple

Second-stage regression using fitted values from first stage.

weak_instrument_test(model::LPIVModel{T}; threshold::T=T(10.0)) -> NamedTuple

Test for weak instruments using Stock-Yogo rule of thumb (F > 10).

bspline_basis(horizons::AbstractVector{Int}, degree::Int, n_interior_knots::Int;
              T::Type{<:AbstractFloat}=Float64) -> BSplineBasis{T}

Construct B-spline basis matrix for given horizons.

bspline_basis_value(x::T, i::Int, degree::Int, knots::Vector{T}) where T

Evaluate B-spline basis function using Cox-de Boor recursion.

compare_smooth_lp(Y::AbstractMatrix{T}, shock_var::Int, horizon::Int;
                  lambda::T=T(1.0), kwargs...) -> NamedTuple

Compare standard LP and smooth LP.

cross_validate_lambda(Y::AbstractMatrix{T}, shock_var::Int, horizon::Int;
                      lambda_grid::Vector{T}=T.(10.0 .^ (-4:0.5:2)),
                      k_folds::Int=5, kwargs...) -> T

Select optimal λ via k-fold cross-validation.

estimate_smooth_lp(Y::AbstractMatrix{T}, shock_var::Int, horizon::Int;
                   degree::Int=3, n_knots::Int=4, lambda::T=T(0.0),
                   lags::Int=4, response_vars::Vector{Int}=collect(1:size(Y,2)),
                   cov_type::Symbol=:newey_west, bandwidth::Int=0) -> SmoothLPModel{T}

Estimate Smooth LP with B-spline parameterization (Barnichon & Brownlees 2019).

roughness_penalty_matrix(basis::BSplineBasis{T}) -> Matrix{T}

Compute roughness penalty matrix R for B-splines (second derivative penalty).

smooth_lp_irf(model::SmoothLPModel{T}; conf_level::Real=0.95) -> LPImpulseResponse{T}

Extract smoothed IRF from SmoothLPModel.

estimate_state_lp(Y::AbstractMatrix{T}, shock_var::Int, state_var::AbstractVector{T},
                  horizon::Int; gamma::Union{T,Symbol}=:estimate, ...) -> StateLPModel{T}

Estimate state-dependent LP (Auerbach & Gorodnichenko 2013).

estimate_transition_params(state_var::AbstractVector{T}, Y::AbstractMatrix{T},
                           shock_var::Int; method::Symbol=:nlls, ...) -> NamedTuple

Estimate smooth transition parameters (γ, c).

exponential_transition(z::AbstractVector{T}, gamma::T, c::T) -> Vector{T}

Exponential (symmetric) transition: F(z) = 1 - exp(-γ(z - c)²)

indicator_transition(z::AbstractVector{T}, c::T) -> Vector{T}

Sharp indicator transition: F(z) = 1 if z ≥ c, else 0

logistic_transition(z::AbstractVector{T}, gamma::T, c::T) -> Vector{T}

Logistic transition function: F(z) = exp(-γ(z - c)) / (1 + exp(-γ(z - c)))

state_irf(model::StateLPModel{T}; regime::Symbol=:both, conf_level::Real=0.95) -> NamedTuple

Extract state-dependent IRFs.

test_regime_difference(model::StateLPModel{T}; h::Union{Int,Nothing}=nothing) -> NamedTuple

Test whether IRFs differ across regimes.

doubly_robust_lp(Y::AbstractMatrix{T}, treatment::AbstractVector{Bool},
                 covariates::AbstractMatrix{T}, horizon::Int; ...) -> PropensityLPModel{T}

Doubly robust LP estimator combining IPW and regression adjustment.

estimate_propensity_lp(Y::AbstractMatrix{T}, treatment::AbstractVector{Bool},
                       covariates::AbstractMatrix{T}, horizon::Int; ...) -> PropensityLPModel{T}

Estimate LP with Inverse Propensity Weighting (Angrist et al. 2018).

estimate_propensity_score(treatment::AbstractVector{Bool}, X::AbstractMatrix{T};
                          method::Symbol=:logit) -> Vector{T}

Estimate propensity scores P(D=1|X) via logit or probit.

inverse_propensity_weights(treatment::AbstractVector{Bool}, propensity::AbstractVector{T};
                           trimming::Tuple{T,T}=(T(0.01), T(0.99)), normalize::Bool=true) -> Vector{T}

Compute IPW weights with optional trimming.

propensity_diagnostics(model::PropensityLPModel{T}) -> NamedTuple

Propensity score diagnostics (overlap, balance).

propensity_irf(model::PropensityLPModel{T}; conf_level::Real=0.95) -> LPImpulseResponse{T}

Extract treatment effect (ATE) IRF from PropensityLPModel.

Build the control vector from LP model's last observations.

Compute forecast CIs from SEs using normal quantiles.

Bootstrap CIs for LP forecasts via residual resampling.

forecast(lp::LPModel{T}, shock_path::AbstractVector{<:Real};
         ci_method=:analytical, conf_level=0.95, n_boot=500) -> LPForecast{T}

Compute direct multi-step LP forecasts given a shock trajectory.

For each horizon h=1,...,H, the forecast uses the LP regression coefficients: ŷ{T+h} = αh + βh·shockh + Γh·controlsT

where controls_T are the last lags observations of Y.

Arguments

• lp: Estimated LP model
• shock_path: Vector of length H with assumed future shock values

Keyword Arguments

• ci_method: :analytical (default), :bootstrap, or :none
• conf_level: Confidence level (default: 0.95)
• n_boot: Number of bootstrap replications (default: 500)

Returns

LPForecast{T} with point forecasts, CIs, and standard errors.

forecast(slp::StructuralLP{T}, shock_idx::Int, shock_path::AbstractVector{<:Real};
         ci_method=:analytical, conf_level=0.95, n_boot=500) -> LPForecast{T}

Compute direct multi-step forecast from a structural LP model using a specific orthogonalized shock.

Arguments

• slp: Structural LP model
• shock_idx: Index of the structural shock to use (1:n)
• shock_path: Vector of length H with assumed shock values

Keyword Arguments

• ci_method: :analytical (default), :bootstrap, or :none
• conf_level: Confidence level (default: 0.95)
• n_boot: Number of bootstrap replications (default: 500)

Dispatch to R², LP-A, or LP-B estimator at horizon h.

Bootstrap bias correction and CIs for LP-FEVD.

1. Fit bivariate VAR(L) on (z, y) with HQIC lag selection
2. Compute 'true' FEVD from VAR (theoretical benchmark)
3. Simulate B samples from VAR, compute LP-FEVD for each
4. Bias = mean(bootstrap FEVD) - true FEVD
5. Bias-corrected = raw - bias
6. CIs from centered bootstrap distribution (Kilian 1998)

_lp_fevd_lpa_h(lp_model, shock, resp_idx, h) -> T

LP-A FEVD: ŝh = (Σ{i=0}^{h} β̂₀^{i,LP}² σ̂z²) / Var(f̂{t+h|t-1}). Uses IRF coefficients directly — no R² regression needed.

_lp_fevd_lpb_h(f_hat, shock, lp_model, resp_idx, t_start, h) -> T

LP-B FEVD: ŝh = numerator / (numerator + Var(ṽ{t+h|t-1})), where ṽ is the residual from the R²-regression (Eq. 6).

_lp_fevd_r2_h(f_hat, shock, t_start, h) -> T

R² FEVD at horizon h: regress forecast errors f̂ on shock leads [z{t+h},...,zt].

Regress fhat on [1, z{t+h}, z{t+h-1}, ..., zt] and return R².

Regress f_hat on shock leads and return Var(residuals).

Compute R²-based FEVD for scalar (z, y) pair at horizon h. Used in bootstrap bias correction.

Simulate T_sim observations from a VAR model with burn-in.

Compute theoretical FEVD of variable 2 w.r.t. shock 1 from bivariate VAR.

fevd(slp::StructuralLP{T}, horizon::Int; kwargs...) -> LPFEVD{T}

Compute LP-based FEVD for structural LP results. Dispatches to lp_fevd.

See lp_fevd for full documentation.

lp_fevd(slp::StructuralLP{T}, horizon::Int; kwargs...) -> LPFEVD{T}

Compute LP-based FEVD using the R²-based estimator of Gorodnichenko & Lee (2019).

At each horizon h, the share of variable i's forecast error variance due to shock j is estimated by regressing LP forecast error residuals on structural shock leads z{t+h}, z{t+h-1}, ..., z_t and computing R².

Arguments

• slp: Structural LP result from structural_lp()
• horizon: Maximum FEVD horizon (capped at IRF horizon)

Keyword Arguments

• method: Estimator (:r2, :lpa, :lpb). Default: :r2
• bias_correct: Apply VAR-based bootstrap bias correction. Default: true
• n_boot: Number of bootstrap replications. Default: 500
• conf_level: Confidence level for CIs. Default: 0.95
• var_lags: VAR lag order for bias correction (default: HQIC-selected)

Returns

LPFEVD{T} with raw proportions, bias-corrected values, SEs, and CIs.

Reference

Gorodnichenko, Y. & Lee, B. (2019). "Forecast Error Variance Decompositions with Local Projections." JBES, 38(4), 921–933.

_estimate_restricted_em(X, r, blocks, standardize; max_iter=500, tol=1e-6) -> FactorModel

Estimate block-restricted factor model via EM algorithm.

Each factor loads only on its assigned block of variables. The restriction mask R enforces zero loadings outside each block.

Algorithm

1. Validate block structure (count, overlap, range, minimum size)
2. Build N × r restriction mask R from block assignments
3. Initialize Λ via block-wise PCA (first eigenvector per block)
4. EM iteration:
   • E-step: F = X Λ (Λ'Λ)⁻¹ (posterior mean of factors given loadings)
   • M-step: Λ_new = (F'F)⁻¹ F'X, masked by R (zero out restricted entries)
5. Convergence when max absolute change in Λ < tol

estimate_factors(X, r; standardize=true, blocks=nothing) -> FactorModel

Estimate static factor model Xt = Λ Ft + e_t via Principal Component Analysis, or via block-restricted EM when blocks is provided.

Arguments

• X: Data matrix (T × N), observations × variables
• r: Number of factors to extract

Keyword Arguments

• standardize::Bool=true: Standardize data before estimation
• blocks::Union{Nothing, Dict{Symbol, Vector{Int}}}=nothing: Block restriction map. When provided, each key is a block (factor) name and the value is a vector of variable indices that load on that factor. Variables not in a block have zero loadings on that factor. The number of blocks must equal r.

Returns

FactorModel containing factors, loadings, eigenvalues, and explained variance. When blocks is provided, block_names is set on the returned model.

Example

X = randn(200, 50)  # 200 observations, 50 variables
fm = estimate_factors(X, 3)  # Extract 3 factors via PCA
r2(fm)  # R² for each variable

# Block-restricted estimation
blocks = Dict(:real => [1,2,3,4,5], :nominal => [6,7,8,9,10])
fm_restricted = estimate_factors(X, 2; blocks=blocks)

forecast(model::FactorModel, h; p=1, ci_method=:theoretical, conf_level=0.95, n_boot=1000)

Forecast factors and observables h steps ahead from a static factor model.

Internally fits a VAR(p) on the extracted factors, then uses the VAR dynamics to produce multi-step forecasts and confidence intervals.

Arguments

• model: Estimated static factor model
• h: Forecast horizon

Keyword Arguments

• p::Int=1: VAR lag order for factor dynamics
• ci_method::Symbol=:theoretical: CI method — :none, :theoretical, or :bootstrap
• conf_level::Real=0.95: Confidence level for intervals
• n_boot::Int=1000: Number of bootstrap replications (if ci_method=:bootstrap)

Returns

FactorForecast with factor and observable forecasts (and CIs if requested).

ic_criteria(X, max_factors; standardize=true)

Compute Bai-Ng (2002) information criteria IC1, IC2, IC3 for selecting the number of factors.

Arguments

• X: Data matrix (T × N)
• max_factors: Maximum number of factors to consider

Returns

Named tuple with IC values and optimal factor counts:

• IC1, IC2, IC3: Information criteria vectors
• r_IC1, r_IC2, r_IC3: Optimal factor counts

Example

result = ic_criteria(X, 10)
println("Optimal factors: IC1=", result.r_IC1, ", IC2=", result.r_IC2, ", IC3=", result.r_IC3)

scree_plot_data(m::FactorModel)

Return data for scree plot: factor indices, explained variance, cumulative variance.

Example

data = scree_plot_data(fm)
# Plot: data.factors vs data.explained_variance

Degrees of freedom.

Number of observations.

Predicted values: F * Λ'.

R² for each variable.

Residuals: X - predicted.

Compute Gaussian log-likelihood given factors and parameters.

M-step: update loadings, VAR coefficients, and covariances.

EM algorithm for maximum likelihood estimation.

Two-step estimation: PCA for factors, then VAR on extracted factors.

companion_matrix_factors(model::DynamicFactorModel)

Construct companion matrix for factor VAR dynamics.

The companion form converts VAR(p) to VAR(1): [Ft; F{t-1}; ...] = C [F{t-1}; F{t-2}; ...] + [η_t; 0; ...]

estimate_dynamic_factors(X, r, p; method=:twostep, standardize=true, max_iter=100, tol=1e-6)

Estimate dynamic factor model with VAR(p) factor dynamics.

Arguments

• X: Data matrix (T × N)
• r: Number of factors
• p: Number of lags in factor VAR

Keyword Arguments

• method::Symbol=:twostep: Estimation method (:twostep or :em)
• standardize::Bool=true: Standardize data
• max_iter::Int=100: Maximum EM iterations (if method=:em)
• tol::Float64=1e-6: Convergence tolerance (if method=:em)
• diagonal_idio::Bool=true: Assume diagonal idiosyncratic covariance

Returns

DynamicFactorModel with factors, loadings, VAR coefficients, and diagnostics.

Example

dfm = estimate_dynamic_factors(X, 3, 2)  # 3 factors, VAR(2)
forecast(dfm, 12)  # 12-step ahead forecast

forecast(model::DynamicFactorModel, h; ci_method=:theoretical, conf_level=0.95, n_boot=1000, ci=false, ci_level=0.95)

Forecast factors and observables h steps ahead.

Arguments

• model: Estimated dynamic factor model
• h: Forecast horizon

Keyword Arguments

• ci_method::Symbol=:theoretical: CI method — :none, :theoretical, :bootstrap, or :simulation
• conf_level::Real=0.95: Confidence level for intervals
• n_boot::Int=1000: Bootstrap replications (for :bootstrap and :simulation)
• ci::Bool=false: Legacy keyword — ci=true maps to ci_method=:simulation
• ci_level::Real=0.95: Legacy keyword — maps to conf_level

Returns

FactorForecast with factor and observable forecasts (and CIs if requested).

Example

fc = forecast(dfm, 12; ci_method=:theoretical)
fc.observables       # h×N matrix of forecasts
fc.observables_lower # h×N lower CI bounds

ic_criteria_dynamic(X, max_r, max_p; standardize=true, method=:twostep)

Select (r, p) via AIC/BIC grid search over factor and lag combinations.

Returns

Named tuple with AIC/BIC matrices and optimal (r, p) combinations.

is_stationary(model::DynamicFactorModel) -> Bool

Check if factor dynamics are stationary (max |eigenvalue| < 1).

Akaike Information Criterion.

Bayesian Information Criterion.

Degrees of freedom: Nr + r²p + r(r+1)/2 + N.

Log-likelihood of the fitted model.

Number of observations.

Predicted values: F * Λ'.

R² for each variable.

Residuals: X - predicted.

Compute spectral density of common component from loadings and eigenvalues.

Compute kernel weights for spectral smoothing.

Compute variance explained by first q factors (averaged across frequencies).

Estimate spectral density matrix with kernel smoothing.

Extract time-domain factors via frequency-domain projection.

AR(1) forecast for each factor series.

Reconstruct common component in time domain via inverse FFT.

Automatic bandwidth selection: T^(1/3).

Eigendecomposition of spectral density at each frequency.

common_variance_share(model::GeneralizedDynamicFactorModel) -> Vector

Fraction of each variable's variance explained by the common component.

estimate_gdfm(X, q; standardize=true, bandwidth=0, kernel=:bartlett, r=0) -> GeneralizedDynamicFactorModel

Estimate Generalized Dynamic Factor Model using spectral methods.

Arguments

• X: Data matrix (T × N)
• q: Number of dynamic factors

Keyword Arguments

• standardize::Bool=true: Standardize data
• bandwidth::Int=0: Kernel bandwidth (0 = automatic selection)
• kernel::Symbol=:bartlett: Kernel for spectral smoothing (:bartlett, :parzen, :tukey)
• r::Int=0: Number of static factors (0 = same as q)

Returns

GeneralizedDynamicFactorModel with common/idiosyncratic components and spectral loadings.

Example

gdfm = estimate_gdfm(X, 3)
common_variance_share(gdfm)  # Fraction of variance explained by common component

forecast(model::GeneralizedDynamicFactorModel, h; method=:ar, ci_method=:theoretical, conf_level=0.95, n_boot=1000)

Forecast h steps ahead using AR extrapolation of factors.

Arguments

• model: Estimated GDFM
• h: Forecast horizon

Keyword Arguments

• method::Symbol=:ar: Forecasting method (currently only :ar supported)
• ci_method::Symbol=:theoretical: CI method — :none, :theoretical, or :bootstrap
• conf_level::Real=0.95: Confidence level for intervals
• n_boot::Int=1000: Bootstrap replications (for :bootstrap)

Returns

FactorForecast with factor and observable forecasts (and CIs if requested).

ic_criteria_gdfm(X, max_q; standardize=true, bandwidth=0, kernel=:bartlett)

Information criteria for selecting number of dynamic factors.

Uses eigenvalue ratio test and cumulative variance threshold.

Returns

Named tuple with:

• eigenvalue_ratios: Ratios of consecutive eigenvalues
• cumulative_variance: Cumulative variance explained
• q_ratio: Optimal q from eigenvalue ratio
• q_variance: Optimal q from 90% variance threshold

spectral_eigenvalue_plot_data(model::GeneralizedDynamicFactorModel)

Return data for plotting eigenvalues across frequencies.

Degrees of freedom.

Number of observations.

Predicted values (common component).

R² for each variable.

Residuals (idiosyncratic component).

estimate_pvar(d::PanelData{T}, p::Int; kwargs...) -> PVARModel{T}

Estimate Panel VAR(p) via GMM.

Arguments

• d::PanelData{T} — balanced panel data (use xtset to construct)
• p::Int — number of lags

Keyword Arguments

• dependent_vars::Union{Vector{String},Nothing}=nothing — endogenous variable names (default: all)
• predet_vars::Vector{String}=String[] — predetermined variable names
• exog_vars::Vector{String}=String[] — strictly exogenous variable names
• transformation::Symbol=:fd — :fd (first-difference) or :fod (forward orthogonal deviations)
• steps::Symbol=:twostep — :onestep, :twostep, or :mstep (iterated)
• system_instruments::Bool=false — if true, use System GMM (Blundell-Bond)
• system_constant::Bool=true — include constant in level equation (System GMM)
• min_lag_endo::Int=2 — minimum instrument lag for endogenous variables
• max_lag_endo::Int=99 — maximum instrument lag (99 = all available)
• collapse::Bool=false — collapse instruments to limit proliferation
• pca_instruments::Bool=false — apply PCA reduction to instruments
• pca_max_components::Int=0 — max PCA components (0 = auto)
• max_iter::Int=100 — max iterations for iterated GMM

Returns

• PVARModel{T} with coefficient estimates, robust standard errors, and GMM internals

References

• Arellano, M. & Bond, S. (1991). Review of Economic Studies 58(2), 277-297.
• Blundell, R. & Bond, S. (1998). Journal of Econometrics 87(1), 115-143.
• Windmeijer, F. (2005). Journal of Econometrics 126(1), 25-51.

Examples

pd = xtset(df, :id, :time)
model = estimate_pvar(pd, 2; steps=:twostep)
model = estimate_pvar(pd, 1; system_instruments=true, steps=:twostep)

estimate_pvar_feols(d::PanelData{T}, p::Int; kwargs...) -> PVARModel{T}

Estimate Panel VAR(p) via Fixed Effects OLS (within estimator).

Simpler than GMM but subject to Nickell (1981) bias when T is small relative to N. Uses within-group demeaning to remove fixed effects, then pooled OLS with cluster-robust standard errors.

Arguments

• d::PanelData{T} — panel data
• p::Int — number of lags

Keyword Arguments

• dependent_vars::Union{Vector{String},Nothing}=nothing — endogenous variable names
• predet_vars::Vector{String}=String[] — predetermined variable names
• exog_vars::Vector{String}=String[] — strictly exogenous variable names

Examples

pd = xtset(df, :id, :time)
model = estimate_pvar_feols(pd, 2)

pvar_oirf(model::PVARModel{T}, H::Int) -> Array{T, 3}

Orthogonalized impulse response functions via Cholesky decomposition.

Ψh = Φh P where P = chol(Σ) is the lower Cholesky factor.

Returns H+1 × m × m array: oirf[h+1, response, shock].

Arguments

• model::PVARModel — estimated PVAR
• H::Int — maximum horizon

Examples

oirf = pvar_oirf(model, 10)
oirf[1, :, :]  # impact response (h=0)

pvar_girf(model::PVARModel{T}, H::Int) -> Array{T, 3}

Generalized impulse response functions (Pesaran & Shin, 1998).

GIRFh(j) = Φh Σ ej / √σjj

where ej is the j-th unit vector and σjj = Σ[j,j].

Returns H+1 × m × m array: girf[h+1, response, shock].

pvar_fevd(model::PVARModel{T}, H::Int) -> Array{T, 3}

Forecast error variance decomposition based on OIRF.

Ω[l, k, h] = Σ{j=0}^{h} (Ψj[l,k])² / MSE_h[l,l]

where MSEh = Σ{j=0}^{h} Φj Σ Φj'.

Returns H+1 × m × m array: fevd[h+1, variable, shock]. Each row sums to 1.

Examples

fv = pvar_fevd(model, 10)
sum(fv[11, 1, :])  # ≈ 1.0 (all shocks for var 1 at h=10)

pvar_stability(model::PVARModel{T}) -> PVARStability{T}

Check stability of the PVAR by computing eigenvalues of the companion matrix. The system is stable if all eigenvalue moduli are strictly less than 1.

Examples

stab = pvar_stability(model)
stab.is_stable  # true if all |λ| < 1
stab.moduli     # sorted eigenvalue moduli

pvar_bootstrap_irf(model::PVARModel{T}, H::Int;
                    irf_type::Symbol=:oirf, n_draws::Int=500,
                    ci::Real=0.95, rng::AbstractRNG=Random.default_rng()) -> NamedTuple

Group-level block bootstrap for PVAR impulse response confidence intervals.

Resamples N groups with replacement → re-estimates PVAR → computes IRF → collects pointwise quantiles.

Arguments

• model::PVARModel — estimated PVAR
• H::Int — maximum IRF horizon

Keywords

• irf_type::Symbol=:oirf — :oirf or :girf
• n_draws::Int=500 — number of bootstrap replications
• ci::Real=0.95 — confidence level
• rng::AbstractRNG — random number generator

Returns

NamedTuple with:

• irf::Array{T,3} — point estimate (H+1 × m × m)
• lower::Array{T,3} — lower CI bound
• upper::Array{T,3} — upper CI bound
• draws::Array{T,4} — all bootstrap draws (n_draws × H+1 × m × m)

pvar_hansen_j(model::PVARModel{T}) -> PVARTestResult{T}

Hansen (1982) J-test for overidentifying restrictions.

J = (Σi Zi' ei)' W (Σi Zi' ei) ~ χ²(q - k)

where q = number of instruments, k = number of parameters per equation.

H0: All moment conditions are valid. H1: Some moment conditions are invalid.

Examples

j = pvar_hansen_j(model)
j.pvalue > 0.05  # fail to reject → instruments valid

pvar_mmsc(model::PVARModel{T}; hq_criterion::Real=2.1) -> NamedTuple

Andrews-Lu (2001) Model and Moment Selection Criteria based on Hansen J-statistic.

MMSCBIC = J - (c - b) × log(n) MMSCAIC = J - (c - b) × 2 MMSC_HQIC = J - Q(c - b) × log(log(n))

Lower values are preferred.

Returns

NamedTuple (bic, aic, hqic) of MMSC values.

Examples

mmsc = pvar_mmsc(model)
mmsc.bic   # MMSC-BIC value

pvar_lag_selection(d::PanelData{T}, max_p::Int; kwargs...) -> NamedTuple

Select optimal PVAR lag order by estimating models for p = 1, ..., max_p and comparing Andrews-Lu MMSC criteria.

Returns

NamedTuple with:

• table::Matrix — comparison table (max_p × 4: p, BIC, AIC, HQIC)
• best_bic::Int, best_aic::Int, best_hqic::Int — optimal lag orders
• models::Vector{PVARModel} — estimated models

Examples

sel = pvar_lag_selection(pd, 4)
sel.best_bic  # optimal lag by BIC

linear_gmm_solve(S_ZX::Matrix{T}, S_Zy::AbstractVecOrMat{T},
                  W::Matrix{T}) -> Vector{T}

Solve linear GMM analytically: E[Z'(y - Xβ)] = 0.

Given pre-aggregated cross-products SZX = Σi Zi'Xi and SZy = Σi Zi'yi:

β = (S_ZX' W S_ZX)⁻¹ S_ZX' W S_Zy

Arguments

• S_ZX — q × k aggregated instrument-regressor cross-product
• S_Zy — q × 1 (or q-vector) aggregated instrument-response cross-product
• W — q × q weighting matrix

Returns

• Parameter vector β (length k)

gmm_sandwich_vcov(S_ZX::Matrix{T}, W::Matrix{T},
                   D_e::Matrix{T}) -> Matrix{T}

Robust sandwich covariance for one-step GMM:

V = (S_ZX' W S_ZX)⁻¹ S_ZX' W D_e W S_ZX (S_ZX' W S_ZX)⁻¹

where De = Σi (Zi ei)(Zi ei)' is the robust moment covariance.

Arguments

• S_ZX — q × k aggregated instrument-regressor cross-product
• W — q × q weighting matrix
• D_e — q × q robust moment covariance

Returns

• k × k covariance matrix

andrews_lu_mmsc(j_stat::T, n_instruments::Int, n_params::Int,
                 n_obs::Int; hq_criterion::Real=2.1) -> NamedTuple

Andrews-Lu (2001) Model and Moment Selection Criteria based on Hansen J-statistic.

MMSC_BIC  = J - (c - b) × log(n)
MMSC_AIC  = J - (c - b) × 2
MMSC_HQIC = J - Q(c - b) × log(log(n))

where c = moment conditions, b = parameters, n = observations. Lower values indicate better model specification.

Arguments

• j_stat — Hansen J-test statistic
• n_instruments — number of moment conditions (c)
• n_params — number of parameters (b)
• n_obs — number of observations (n)
• hq_criterion — penalty parameter Q for HQIC (default 2.1)

Returns

NamedTuple (bic, aic, hqic).

References

• Andrews, D. & Lu, B. (2001). Journal of Econometrics 101(1), 123-164.

estimate_did(pd::PanelData{T}, outcome::Union{String,Symbol},
             treatment::Union{String,Symbol};
             method::Symbol=:twfe,
             leads::Int=0, horizon::Int=5,
             covariates::Vector{String}=String[],
             control_group::Symbol=:never_treated,
             cluster::Symbol=:unit,
             conf_level::Real=0.95) where {T}

Estimate a Difference-in-Differences model.

Arguments

• pd: Panel data with outcome and treatment timing columns
• outcome: Name of outcome variable
• treatment: Name of treatment timing variable (contains period of first treatment; 0 or NaN for never-treated units)
• method: Estimation method
  • :twfe – Two-Way Fixed Effects event-study regression (default)
  • :callaway_santanna – Callaway & Sant'Anna (2021) group-time ATT
  • :sun_abraham – Sun & Abraham (2021) interaction-weighted estimator
  • :bjs – Borusyak, Jaravel & Spiess (2024) imputation estimator
  • :did_multiplegt – de Chaisemartin & D'Haultfoeuille (2020)

Keyword Arguments

• leads: Number of pre-treatment periods to estimate (default: 0)
• horizon: Post-treatment horizon (default: 5)
• covariates: Additional control variable names
• control_group: :never_treated (default) or :not_yet_treated
• cluster: SE clustering: :unit (default), :time, :twoway
• conf_level: Confidence level (default: 0.95)
• n_boot: Number of bootstrap replications for :did_multiplegt (default: 200)
• base_period: :varying (default, R did default) or :universal. For :callaway_santanna only. :varying uses adjacent-period comparisons for pre-treatment cells; :universal always uses g-1 as base period

Returns

DIDResult{T} – unified result type for all methods.

Examples

# TWFE event study
did = estimate_did(pd, :gdp, :reform_year; method=:twfe, leads=3, horizon=5)
report(did)

# Callaway-Sant'Anna
did_cs = estimate_did(pd, :gdp, :reform_year; method=:callaway_santanna, leads=3, horizon=5)
plot_result(did_cs)

References

• Callaway, B. & Sant'Anna, P. H. C. (2021). JoE 225(2), 200-230.

estimate_event_study_lp(pd::PanelData{T}, outcome::Union{String,Symbol},
                        treatment::Union{String,Symbol}, H::Int;
                        leads::Int=3, lags::Int=4,
                        covariates::Vector{String}=String[],
                        cluster::Symbol=:unit,
                        conf_level::Real=0.95) where {T}

Estimate event study impulse responses using local projections with panel FE.

Uses the switching indicator ΔD_{it} (= 1 only at the treatment onset period) and time-only FE. Already-treated observations are excluded from the sample.

Arguments

• pd: Panel data
• outcome: Outcome variable name
• treatment: Treatment variable (binary 0/1 or timing column)
• H: Maximum post-treatment horizon
• leads: Number of pre-treatment leads to estimate (for pre-trends)
• lags: Number of outcome lags as controls
• covariates: Additional control variable names
• cluster: Clustering level (:unit, :time, :twoway)
• conf_level: Confidence level (default: 0.95)

Returns

EventStudyLP{T} with dynamic treatment effect coefficients.

References

• Jordà, Ò. (2005). AER 95(1), 161-182.
• Acemoglu, D. et al. (2019). JPE 127(1), 47-100.

estimate_lp_did(pd::PanelData{T}, outcome, treatment, H; kwargs...) -> LPDiDResult{T}

LP-DiD estimator with full feature parity to the Stata lpdid package v1.0.2.

Uses the switching indicator (ΔD_{it}) as treatment, clean control samples (CCS), and time-only fixed effects (long differencing absorbs unit FE).

Arguments

• pd — panel data
• outcome — outcome variable (name or symbol)
• treatment — treatment variable (binary 0/1 or timing column)
• H — maximum horizon

Keyword Arguments

• pre_window::Int=3 — pre-treatment event-time window
• post_window::Int=H — post-treatment event-time window
• ylags::Int=0 — number of outcome lags as controls
• dylags::Int=0 — number of ΔY lags as controls
• covariates::Vector{String}=String[] — additional covariates
• nonabsorbing::Union{Nothing,Int}=nothing — stabilization window L for non-absorbing treatment
• notyet::Bool=false — restrict controls to not-yet-treated
• nevertreated::Bool=false — restrict controls to never-treated
• firsttreat::Bool=false — use only first treatment event per unit
• oneoff::Bool=false — one-off treatment (requires nonabsorbing)
• pmd::Union{Nothing,Symbol,Int}=nothing — pre-mean differencing (:max or integer k)
• reweight::Bool=false — IPW reweighting for equally weighted ATE
• nocomp::Bool=false — restrict to observations in CCS at all horizons
• cluster::Symbol=:unit — clustering level (:unit, :time, :twoway)
• conf_level::Real=0.95 — confidence level
• only_pooled::Bool=false — skip event study, compute pooled only
• only_event::Bool=false — skip pooled estimates
• post_pooled::Union{Nothing,Tuple{Int,Int}}=nothing — pooled post-treatment window (start, end)
• pre_pooled::Union{Nothing,Tuple{Int,Int}}=nothing — pooled pre-treatment window (start, end)

References

• Dube, A., Girardi, D., Jordà, Ò. & Taylor, A.M. (2025). JAE.
• Acemoglu, D. et al. (2019). JPE 127(1), 47-100.

bacon_decomposition(pd::PanelData{T}, outcome::Union{String,Symbol},
                    treatment::Union{String,Symbol}) where {T}

Decompose the TWFE DiD estimator into a weighted average of all 2x2 DiD comparisons.

Each weight reflects the subsample size and variance of treatment within the comparison.

Returns

BaconDecomposition{T} with estimates, weights, and comparison types.

Reference

Goodman-Bacon, A. (2021). "Difference-in-Differences with Variation in Treatment Timing." JoE 225(2), 254-277.

pretrend_test(result::DIDResult{T}) where {T}
pretrend_test(result::EventStudyLP{T}) where {T}

Joint F-test for parallel trends: H0: beta{-K} = beta{-K+1} = ... = beta_{-1} = 0.

Tests whether pre-treatment event-time coefficients are jointly zero. High p-value -> no evidence against parallel trends.

Returns

PretrendTestResult{T} with F-statistic and p-value.

negative_weight_check(pd::PanelData{T}, treatment::Union{String,Symbol}) where {T}

de Chaisemartin & D'Haultfoeuille (2020) negative weight diagnostic.

Checks whether the TWFE estimator assigns negative weights to some group-time cells. Negative weights can cause the overall DiD estimate to have the opposite sign of all underlying ATT(g,t).

Returns

NegativeWeightResult{T}.

Reference

de Chaisemartin, C. & D'Haultfoeuille, X. (2020). AER 110(9), 2964-2996.

honest_did(result::DIDResult{T}; Mbar::Real=1.0, conf_level::Real=0.95) where {T}

HonestDiD sensitivity analysis for a Difference-in-Differences result.

Constructs robust confidence intervals that account for bounded violations of the parallel trends assumption under the relative magnitudes restriction (Rambachan & Roth 2023).

Arguments

• result::DIDResult{T} — DiD estimation result (from estimate_did or callaway_santanna)
• Mbar::Real=1.0 — maximum violation bound (relative magnitudes parameter)
• conf_level::Real=0.95 — confidence level for the robust CIs

Returns

HonestDiDResult{T} with robust CIs, breakdown value, and post-treatment estimates.

Example

did = estimate_did(pd, :y, :treat)
h = honest_did(did; Mbar=1.0)
h.breakdown_value  # smallest Mbar that overturns significance

honest_did(result::EventStudyLP{T}; Mbar::Real=1.0, conf_level::Real=0.95) where {T}

HonestDiD sensitivity analysis for an Event Study LP result.

Constructs robust confidence intervals that account for bounded violations of the parallel trends assumption under the relative magnitudes restriction (Rambachan & Roth 2023).

Arguments

• result::EventStudyLP{T} — event study LP result (from event_study_lp or lp_did)
• Mbar::Real=1.0 — maximum violation bound (relative magnitudes parameter)
• conf_level::Real=0.95 — confidence level for the robust CIs

Returns

HonestDiDResult{T} with robust CIs, breakdown value, and post-treatment estimates.

Example

eslp = event_study_lp(pd, :y, :treat, 5; leads=3, lags=2)
h = honest_did(eslp; Mbar=0.5)
h.robust_ci_lower  # lower bounds accounting for trend violations

andrews_lu_mmsc(j_stat::T, n_instruments::Int, n_params::Int,
                 n_obs::Int; hq_criterion::Real=2.1) -> NamedTuple

Andrews-Lu (2001) Model and Moment Selection Criteria based on Hansen J-statistic.

MMSC_BIC  = J - (c - b) × log(n)
MMSC_AIC  = J - (c - b) × 2
MMSC_HQIC = J - Q(c - b) × log(log(n))

where c = moment conditions, b = parameters, n = observations. Lower values indicate better model specification.

Arguments

• j_stat — Hansen J-test statistic
• n_instruments — number of moment conditions (c)
• n_params — number of parameters (b)
• n_obs — number of observations (n)
• hq_criterion — penalty parameter Q for HQIC (default 2.1)

Returns

NamedTuple (bic, aic, hqic).

References

• Andrews, D. & Lu, B. (2001). Journal of Econometrics 101(1), 123-164.

estimate_gmm(moment_fn::Function, theta0::AbstractVector{T}, data;
             weighting::Symbol=:two_step, max_iter::Int=100,
             tol::T=T(1e-8), hac::Bool=true, bandwidth::Int=0,
             bounds::Union{Nothing,ParameterTransform}=nothing) -> GMMModel{T}

Estimate parameters via Generalized Method of Moments.

Minimizes: Q(θ) = g(θ)'W g(θ) where g(θ) = (1/n) Σᵢ gᵢ(θ)

Arguments:

• moment_fn: Function (theta, data) -> Matrix of moment conditions (n × q)
• theta0: Initial parameter guess
• data: Data passed to moment function
• weighting: Weighting method (:identity, :optimal, :two_step, :iterated)
• max_iter: Maximum iterations for optimization and/or iterated GMM
• tol: Convergence tolerance
• hac: Use HAC correction for optimal weighting
• bandwidth: HAC bandwidth (0 = automatic)
• bounds: Parameter bounds via ParameterTransform (default: nothing = unconstrained). When provided, optimization is performed in unconstrained space via bijective transforms, and standard errors are corrected via the delta method.

Returns:

• GMMModel with estimates, covariance, and J-test results

estimate_lp_gmm(Y::AbstractMatrix{T}, shock_var::Int, horizon::Int;
                lags::Int=4, weighting::Symbol=:two_step) -> Vector{GMMModel{T}}

Estimate Local Projection via GMM.

Returns a GMMModel for each horizon.

gmm_objective(theta::AbstractVector{T}, moment_fn::Function, data,
              W::AbstractMatrix{T}) -> T

Compute GMM objective: Q(θ) = g(θ)'W g(θ)

where g(θ) = (1/n) Σᵢ gᵢ(θ,data)

Arguments:

• theta: Parameter vector
• moment_fn: Function (theta, data) -> Matrix of moment conditions (n × q)
• data: Data passed to moment function
• W: Weighting matrix (q × q)

Returns:

• GMM objective value

gmm_sandwich_vcov(S_ZX::Matrix{T}, W::Matrix{T},
                   D_e::Matrix{T}) -> Matrix{T}

Robust sandwich covariance for one-step GMM:

V = (S_ZX' W S_ZX)⁻¹ S_ZX' W D_e W S_ZX (S_ZX' W S_ZX)⁻¹

where De = Σi (Zi ei)(Zi ei)' is the robust moment covariance.

Arguments

• S_ZX — q × k aggregated instrument-regressor cross-product
• W — q × q weighting matrix
• D_e — q × q robust moment covariance

Returns

• k × k covariance matrix

gmm_summary(model::GMMModel{T}) -> NamedTuple

Summary statistics for GMM estimation.

identity_weighting(n_moments::Int, ::Type{T}=Float64) -> Matrix{T}

Identity weighting matrix (one-step GMM).

j_test(model::GMMModel{T}) -> NamedTuple

Hansen's J-test for overidentifying restrictions.

H0: All moment conditions are valid (E[g(θ₀)] = 0) H1: Some moment conditions are violated

Returns:

• J_stat: Test statistic
• p_value: p-value from chi-squared distribution
• df: Degrees of freedom (nmoments - nparams)
• reject_05: Whether to reject at 5% level

linear_gmm_solve(S_ZX::Matrix{T}, S_Zy::AbstractVecOrMat{T},
                  W::Matrix{T}) -> Vector{T}

Solve linear GMM analytically: E[Z'(y - Xβ)] = 0.

Given pre-aggregated cross-products SZX = Σi Zi'Xi and SZy = Σi Zi'yi:

β = (S_ZX' W S_ZX)⁻¹ S_ZX' W S_Zy

Arguments

• S_ZX — q × k aggregated instrument-regressor cross-product
• S_Zy — q × 1 (or q-vector) aggregated instrument-response cross-product
• W — q × q weighting matrix

Returns

• Parameter vector β (length k)

lp_gmm_moments(Y::AbstractMatrix{T}, shock_var::Int, h::Int, theta,
               lags::Int) -> Matrix{T}

Construct moment conditions for LP estimated via GMM.

Moments: E[Zt * ε{t+h}] = 0 where ε{t+h} = y{t+h} - θ' * X_t and Z includes all exogenous variables.

This is useful when you need to impose cross-equation restrictions.

minimize_gmm(moment_fn::Function, theta0::AbstractVector{T}, data,
             W::AbstractMatrix{T}; max_iter::Int=100, tol::T=T(1e-8)) -> NamedTuple

Minimize GMM objective using Optim.jl (LBFGS with NelderMead fallback).

Returns:

• theta: Minimizer
• objective: Final objective value
• converged: Convergence flag
• iterations: Number of iterations

numerical_gradient(f::Function, x::AbstractVector{T}; eps::T=T(1e-7)) -> Matrix{T}

Compute numerical gradient (Jacobian) of function f at point x using central differences.

Arguments:

• f: Function that takes vector x and returns vector (moment conditions)
• x: Point at which to evaluate gradient
• eps: Step size for finite differences

Returns:

• Jacobian matrix (nmoments × nparams)

optimal_weighting_matrix(moment_fn::Function, theta::AbstractVector{T}, data;
                         hac::Bool=true, bandwidth::Int=0) -> Matrix{T}

Compute optimal GMM weighting matrix: W = inv(Var(g)).

For i.i.d. data: W = inv((1/n) Σᵢ gᵢ gᵢ') For time series: Uses HAC estimation with Newey-West kernel.

Arguments:

• moment_fn: Moment function
• theta: Current parameter estimate
• data: Data
• hac: Use HAC correction for serial correlation
• bandwidth: HAC bandwidth (0 = automatic)

Returns:

• Optimal weighting matrix (q × q)

_moment_covariance(data, moments_fn; hac, bandwidth) -> Matrix{T}

Internal: Compute long-run covariance of per-observation moment contributions. Shared by smm_weighting_matrix (inverted) and smm_data_covariance (raw).

autocovariance_moments(data::AbstractMatrix{T}; lags::Int=1) -> Vector{T}

Compute standard DSGE moment vector from data matrix.

Returns: [upper-triangle variance-covariance elements; diagonal autocovariances at each lag]

For k variables, lags L: k*(k+1)/2 + k*L moments total.

Arguments

• data –- T_obs x k data matrix
• lags –- number of autocovariance lags (default: 1)

estimate_smm(simulator_fn, moments_fn, theta0, data;
             sim_ratio=5, burn=100, weighting=:two_step,
             bounds=nothing, hac=true, bandwidth=0,
             max_iter=1000, tol=1e-8,
             rng=Random.default_rng()) -> SMMModel{T}

Estimate parameters via Simulated Method of Moments.

Minimizes Q(theta) = (m_data - m_sim(theta))' W (m_data - m_sim(theta)) where m_data are data moments and m_sim(theta) are simulated moments.

Arguments

• simulator_fn(theta, T_periods, burn; rng) –- simulates T_periods obs after discarding burn
• moments_fn(data) -> Vector{T} –- computes moment vector from any T x k data matrix
• theta0 –- initial parameter guess
• data –- observed data matrix (T x k)

Keywords

• sim_ratio::Int=5 –- tau = simulation periods / data periods
• burn::Int=100 –- burn-in periods for simulator
• weighting::Symbol=:two_step –- :identity or :two_step
• bounds::Union{Nothing,ParameterTransform}=nothing –- parameter bounds
• hac::Bool=true –- HAC for weighting matrix
• bandwidth::Int=0 –- HAC bandwidth (0 = automatic)
• max_iter::Int=1000 –- max optimizer iterations
• tol=1e-8 –- convergence tolerance
• rng –- random number generator (copied inside objective for deterministic simulation)

References

• Ruge-Murcia (2012), Lee & Ingram (1991)

j_test(m::SMMModel{T}) -> NamedTuple

Hansen's J-test for overidentifying restrictions on an SMM model.

H0: All moment conditions are valid (E[g(theta_0)] = 0) H1: Some moment conditions are violated

Returns:

• J_stat: Test statistic
• p_value: p-value from chi-squared distribution
• df: Degrees of freedom (nmoments - nparams)
• reject_05: Whether to reject at 5% level

smm_data_covariance(data::AbstractMatrix{T}, moments_fn::Function;
                      hac::Bool=true, bandwidth::Int=0) -> Matrix{T}

Compute long-run covariance Ω of data moment contributions for sandwich SE formula.

smm_weighting_matrix(data::AbstractMatrix{T}, moments_fn::Function;
                      hac::Bool=true, bandwidth::Int=0) -> Matrix{T}

Compute optimal SMM weighting matrix from data moment contributions. Centers the per-observation moment contributions and applies HAC with Bartlett kernel.

Arguments

• data — T_obs × k data matrix
• moments_fn — function computing moment vector from data
• hac — use HAC correction (default: true)
• bandwidth — HAC bandwidth, 0 = automatic: floor(4*(n/100)^(2/9))

to_unconstrained(pt::ParameterTransform, theta::AbstractVector) -> Vector

Map parameters from constrained (model) space to unconstrained (optimizer) space.

to_constrained(pt::ParameterTransform, phi::AbstractVector) -> Vector

Map parameters from unconstrained (optimizer) space to constrained (model) space.

transform_jacobian(pt::ParameterTransform, phi::AbstractVector) -> Matrix

Diagonal Jacobian ∂θ/∂φ of the inverse transform (unconstrained → constrained). Used for delta method SE correction.

adf_test(y; lags=:aic, max_lags=nothing, regression=:constant) -> ADFResult

Augmented Dickey-Fuller test for unit root.

Tests H₀: y has a unit root (non-stationary) against H₁: y is stationary.

Arguments

• y: Time series vector
• lags: Number of augmenting lags, or :aic/:bic/:hqic for automatic selection
• max_lags: Maximum lags for automatic selection (default: floor(12*(T/100)^0.25))
• regression: Deterministic terms - :none, :constant (default), or :trend

Returns

ADFResult containing test statistic, p-value, critical values, etc.

Example

y = cumsum(randn(200))  # Random walk (has unit root)
result = adf_test(y)
result.pvalue > 0.05  # Should fail to reject H₀

References

• Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. JASA, 74(366), 427-431.
• MacKinnon, J. G. (2010). Critical values for cointegration tests. Queen's Economics Department Working Paper No. 1227.

kpss_test(y; regression=:constant, bandwidth=:auto) -> KPSSResult

Kwiatkowski-Phillips-Schmidt-Shin test for stationarity.

Tests H₀: y is stationary against H₁: y has a unit root.

Arguments

• y: Time series vector
• regression: :constant (level stationarity) or :trend (trend stationarity)
• bandwidth: Bartlett kernel bandwidth, or :auto for Newey-West selection

Returns

KPSSResult containing test statistic, p-value, critical values, etc.

Example

y = randn(200)  # Stationary series
result = kpss_test(y)
result.pvalue > 0.05  # Should fail to reject H₀ (stationarity)

References

• Kwiatkowski, D., Phillips, P. C., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54(1-3), 159-178.

pp_test(y; regression=:constant, bandwidth=:auto) -> PPResult

Phillips-Perron test for unit root with non-parametric correction.

Tests H₀: y has a unit root against H₁: y is stationary.

Arguments

• y: Time series vector
• regression: :none, :constant (default), or :trend
• bandwidth: Newey-West bandwidth, or :auto for automatic selection

Returns

PPResult containing test statistic (Zt), p-value, critical values, etc.

Example

y = cumsum(randn(200))  # Random walk
result = pp_test(y)
result.pvalue > 0.05  # Should fail to reject H₀

References

• Phillips, P. C., & Perron, P. (1988). Testing for a unit root in time series regression. Biometrika, 75(2), 335-346.

AO (additive outlier) model for a single break candidate.

Step 1: Detrend y via OLS on deterministics with break dummies → residuals e. Step 2: ADF regression on e with pulse dummy D(TB)_t = 1{t == tb} and its lags.

Returns (tstatonelag, lags_used).

IO (innovational outlier) model for a single break candidate.

Regression: Δyt = μ + β·t + θ·DUt(tb) [+ δ·DTt(tb)] + γ·y{t-1} + Σδj·Δy{t-j} + et where DUt = 1{t >= tb} and DT_t = max(0, t - tb + 1).

Returns (tstatongamma, lagsused).

za_test(y; regression=:both, trim=0.15, lags=:aic, max_lags=nothing, outlier=:io) -> ZAResult

Zivot-Andrews test for unit root with endogenous structural break.

Tests H₀: y has a unit root without break against H₁: y is stationary with break.

Arguments

• y: Time series vector
• regression: Type of break - :constant (intercept), :trend (slope), or :both
• trim: Trimming fraction for break search (default 0.15)
• lags: Number of augmenting lags, or :aic/:bic for automatic selection
• max_lags: Maximum lags for selection
• outlier: Outlier model - :io (innovational outlier) or :ao (additive outlier)

Returns

ZAResult containing minimum t-statistic, break point, p-value, etc.

Example

# Series with structural break
y = vcat(randn(100), randn(100) .+ 2)
result = za_test(y; regression=:constant)

# Additive outlier model
result_ao = za_test(y; regression=:constant, outlier=:ao)

References

• Zivot, E., & Andrews, D. W. K. (1992). Further evidence on the great crash, the oil-price shock, and the unit-root hypothesis. JBES, 10(3), 251-270.

ngperron_test(y; regression=:constant) -> NgPerronResult

Ng-Perron unit root tests with GLS detrending (MZα, MZt, MSB, MPT).

Tests H₀: y has a unit root against H₁: y is stationary. These tests have better size properties than ADF/PP in small samples.

Arguments

• y: Time series vector
• regression: :constant (default) or :trend

Returns

NgPerronResult containing MZα, MZt, MSB, MPT statistics and critical values.

Example

y = cumsum(randn(100))
result = ngperron_test(y)
# Check if MZt rejects at 5%
result.MZt < result.critical_values[:MZt][5]

References

• Ng, S., & Perron, P. (2001). Lag length selection and the construction of unit root tests with good size and power. Econometrica, 69(6), 1519-1554.

johansen_test(Y, p; deterministic=:constant) -> JohansenResult

Johansen cointegration test for VAR system.

Tests for the number of cointegrating relationships among variables using trace and maximum eigenvalue tests.

Arguments

• Y: Data matrix (T × n)
• p: Number of lags in the VECM representation
• deterministic: Specification for deterministic terms
  • :none - No deterministic terms
  • :constant - Constant in cointegrating relation (default)
  • :trend - Linear trend in levels

Returns

JohansenResult containing trace and max-eigenvalue statistics, cointegrating vectors, adjustment coefficients, and estimated rank.

Example