# BEGE Model Estimation for Structural Shocks

Replication code for estimating time-varying conditional variance and skewness of macroeconomic structural shocks using the BEGE-GARCH framework.

**Reference:** Engstrom, Eric (2025). "Soft Landing or Stagflation? Measuring Macroeconomic Risks in Real Time."

---

## Overview

The **BEGE (Bad Environment, Good Environment)** model represents economic shocks as differences of two gamma distributions with time-varying arrival rates:

```
u(t) ~ Gamma(p(t), θ_p) - Gamma(n(t), θ_n)
```

where:
- **p(t)** = good news arrival rate (time-varying via GARCH)
- **n(t)** = bad news arrival rate (time-varying via GARCH)
- **θ_p, θ_n** = scale parameters

This framework captures time-varying uncertainty in four structural shocks (supply, demand, idiosyncratic unemployment, idiosyncratic inflation) and produces four-quarter-ahead forecast distributions for GDP growth, core inflation, headline inflation, and the unemployment rate.

---

## Data Requirements

**All data are contained in two CSV files:**

### `data_pull_spfdata.csv` (Required)

This CSV file contains quarterly observations of macroeconomic variables and Survey of Professional Forecasters (SPF) forecasts:

**Observed variables:**
- `gdp_obs` - Real GDP growth (four-quarter percent change)
- `core_pce_obs` - Core PCE inflation (four-quarter percent change)
- `pce_obs` - Headline PCE inflation (four-quarter percent change)
- `u_obs` - Unemployment rate (percent)
- `yeffq` - Effective federal funds rate (quarterly average, percent)

**SPF forecasts (four-quarter ahead):**
- `gdp_forecast` - GDP growth forecast
- `pce_forecast` - PCE inflation forecast
- `u_forecast` - Unemployment rate forecast

**Supporting variables:**
- `dates_all` - Date vector (year + quarter/4 format, e.g., 2025.75 = 2025Q3)

**Note:** Missing values are coded as -99 and automatically converted to NaN by the code.

### `data_notes_spfdata.csv` (Documentation)

This file contains detailed documentation for each variable, including:
- Variable descriptions
- Data sources (FRED series codes, SPF variable names)
- Transformations applied
- Time periods covered

**Refer to this file for complete data source documentation and replication details.**

**No other external data sources are required.** The code is self-contained once these two CSV files are provided. NBER recession dates are hardcoded in `RecessionDummy.m`.

---

## Execution Workflow

The code consists of three main programs that **must be run sequentially** in the following order:

### 1. `est_structural_shocks.m` - Extract Structural Shocks

**Purpose:** Identifies four structural shocks from the macro data.

**What it does:**
- Estimates conditional mean dynamics using regression-based forecasts (SPF forecasts + lagged data)
- Identifies structural shocks via GMM on second moments:
  - `rev_sup` - Supply shock
  - `rev_dem` - Demand shock
  - `rev_p` - Idiosyncratic inflation shock
  - `rev_u` - Idiosyncratic unemployment shock
- Creates mapping matrix `M` (structural → reduced-form)
- Generates workspace variables required by subsequent programs

**Outputs:**
- Workspace variables: `rev_sup`, `rev_dem`, `rev_p`, `rev_u`, `M`, `ndate`, `recdum`, `exp_pi_y`
- `Figure1.png` - Conditional means and reduced-form shocks
- `Figure2.png` - Structural shocks

**Configuration:**
- `edate` - End date for sample (line 59, currently set to 2026.25 = 2026Q1)

**To run:**
```matlab
run est_structural_shocks.m
```

---

### 2. `estimate_bege_model.m` - Estimate BEGE Parameters

**Purpose:** Estimates time-varying good/bad news arrival rates for each structural shock.

**What it does:**
- Estimates BEGE-GARCH models for all four structural shocks
- For each shock, tests 4 specifications (static, dynamic p, dynamic n, both dynamic)
- Selects optimal specification using AICc (corrected Akaike Information Criterion)
- Optionally computes bootstrap standard errors
- Calculates conditional moments (variance, skewness) for structural and endogenous variables
- Saves estimated parameters to `outbege.mat` for future use

**Outputs:**
- Workspace variables: `sup_pt`, `sup_nt`, `dem_pt`, `dem_nt`, `p_pt`, `p_nt`, `u_pt`, `u_nt` (arrival rates)
- Workspace variables: `sup_var`, `dem_var`, `u_var`, `p_var` (conditional variances and skewness)
- Workspace variables: `c_var`, `g_var`, `eu_var`, `ep_var` (endogenous variable moments)
- `outbege.mat` - Saved BEGE parameters
- `Figure3.png` - Variance decomposition (good vs bad news)
- `Figure4.png` - Third moments (skewness) of structural shocks
- `Figure8.png` - Variance and skewness of endogenous variables
- `Figure9.png` - Covariance: GDP vs headline inflation

**Configuration (lines 92-100):**
- `re_estimate` - **Set to `true` for first-time use** (line 94)
- `nboot` - Bootstrap iterations: 1000 for publication, 0 to skip (line 97)
- `verbose` - Display progress messages (line 100)
- `lag` - Forecast horizon in quarters (line 104, default = 4)
- `expt` - GARCH exponent (line 106, default = 1)

**To run:**
```matlab
run estimate_bege_model.m
```

---

### 3. `calc_bege_distributions.m` - Generate Forecast Distributions

**Purpose:** Produces bivariate forecast distributions and calculates scenario probabilities.

**What it does:**
- Simulates the quadrivariate BEGE system (100 million draws)
- Creates bivariate contour plots showing forecast distributions:
  - GDP growth vs inflation
  - Unemployment vs inflation
- Calculates scenario probabilities:
  - Soft landing, slowflation, stagflation
  - Quadrant probabilities (high/low growth × high/low inflation)
- Computes distribution of federal funds rate under Taylor (1999) rule
- Generates quantile time series for all macro variables

**Outputs:**
- Workspace probabilities: `prob_softland`, `prob_slowflation`, `prob_stagflation`
- Workspace probabilities: `prob_ghi_pilo`, `prob_ghi_pihi`, `prob_glo_pilo`, `prob_glo_pihi` (GDP-inflation quadrants)
- Workspace probabilities: `prob_uhi_pilo`, `prob_uhi_pihi`, `prob_ulo_pilo`, `prob_ulo_pihi` (unemployment-inflation quadrants)
- Workspace quantiles: `pcts_c`, `pcts_g`, `pcts_u`, `pcts_p` (10th, 25th, 50th, 75th, 90th percentiles)
- Workspace quantiles: `pcts_t99dist` (Taylor rule distribution percentiles)
- `Figure_GDP_Inflation_*.png` - Bivariate distributions (multiple pages, 4 per page)
- `Figure_Unemployment_Inflation_*.png` - Unemployment vs inflation distributions
- `Figure_Taylor99_and_Probabilities.png` - Policy prescriptions and scenario probabilities
- `Figure_Macro_Quantiles.png` - Quantile time series for all variables

**Configuration (lines 58-95):**
- `pistar`, `gstar`, `ustar` - Reference levels for inflation, growth, unemployment (lines 60-62)
- `rstar`, `picoef`, `gapcoef`, `okunc` - Taylor rule coefficients (lines 65-68)
- `nsim` - Number of simulations (line 95, default = 1e8 for production, 1e6 for testing)
- `lookbk` - Number of quarters to plot (line 99, default = 19 quarters = ~5 years)

**To run:**
```matlab
run calc_bege_distributions.m
```

---

## First-Time Setup

For first-time use or full replication, follow this sequence:

```matlab
% Step 1: Extract structural shocks from data
run est_structural_shocks.m

% Step 2: Estimate BEGE parameters
% IMPORTANT: Set re_estimate = true in estimate_bege_model.m (line 94)
% This is required because no .mat files are included in the replication package
run estimate_bege_model.m

% Step 3: Generate forecast distributions
run calc_bege_distributions.m
```

**Configuration for first-time use:**

In `estimate_bege_model.m`, set:
```matlab
re_estimate = true;   % Required for initial estimation (line 94)
nboot = 0;            % Set to 1000 for bootstrap standard errors (line 97)
verbose = true;       % Display progress (line 100)
```

**Why `re_estimate = true`?**
This replication package does not include pre-estimated parameter files (`outbege.mat`). Setting `re_estimate = true` performs full maximum likelihood estimation of all BEGE models, which takes approximately 5-10 minutes (without bootstrap) or 2-4 hours (with `nboot = 1000`).

**Bootstrap standard errors:**
- For testing/replication verification: `nboot = 0` (fast, ~5-10 minutes)
- For publication-quality results: `nboot = 1000` (slow, ~2-4 hours)

---

## BEGE Model Specifications

The code estimates four specifications for each structural shock and selects the best via AICc:

| Spec | Description | Parameters | GARCH Dynamics |
|------|-------------|------------|----------------|
| 1 | Static | `[p0, θ_p, n0, θ_n]` | None (constant pt, nt) |
| 2 | Dynamic p | `[p0, θ_p, φ_p, ρ_p, n0, θ_n]` | `p(t) = p0 + ρ_p*p(t-1) + φ_p*|u(t)|^expt` |
| 3 | Dynamic n | `[p0, θ_p, n0, θ_n, φ_n, ρ_n]` | `n(t) = n0 + ρ_n*n(t-1) + φ_n*|u(t)|^expt` |
| 4 | Both dynamic | `[p0, θ_p, φ_p, ρ_p, n0, θ_n, φ_n, ρ_n]` | Both p(t) and n(t) time-varying |

**AICc selection criterion:**
```
AICc = -2*log(L) + 2k + (2k² + 2k)/(n_eff - k - 1)
```
where `n_eff = T/lag` adjusts for overlapping observations (Hansen & Hodrick, 1980).

---

## Output Summary

### Workspace Variables (Available After Running All Three Programs)

**From est_structural_shocks.m:**
- `rev_sup`, `rev_dem`, `rev_p`, `rev_u` - Structural shocks
- `M` - Mapping matrix (4×4, structural → reduced-form)
- `exp_pi_y` - Conditional mean forecasts
- `ndate`, `yr`, `qt` - Date vectors
- `recdum` - Recession indicator

**From estimate_bege_model.m:**
- `sup_pt`, `sup_nt`, `sup_thp`, `sup_thn` - Supply shock BEGE parameters (and similarly for dem, p, u)
- `sup_var`, `sup_skw` - Supply shock conditional variance and skewness (and similarly for dem, p, u)
- `c_var`, `g_var`, `eu_var`, `ep_var` - Endogenous variable conditional variances
- `c_thr`, `g_thr`, `eu_thr`, `ep_thr` - Endogenous variable third moments
- `cg_cov`, `cu_cov`, `cp_cov`, `gu_cov`, `gp_cov`, `up_cov` - Covariances

**From calc_bege_distributions.m:**
- `prob_softland`, `prob_slowflation`, `prob_stagflation` - Scenario probabilities
- `prob_ghi_pilo`, `prob_ghi_pihi`, `prob_glo_pilo`, `prob_glo_pihi` - GDP-inflation quadrants
- `prob_uhi_pilo`, `prob_uhi_pihi`, `prob_ulo_pilo`, `prob_ulo_pihi` - Unemployment-inflation quadrants
- `pcts_c`, `pcts_g`, `pcts_u`, `pcts_p` - Macro variable quantiles (10th, 25th, 50th, 75th, 90th)
- `pcts_t99dist` - Taylor rule distribution quantiles

### Saved Files

**Parameter files:**
- `outbege.mat` - All estimated BEGE parameters and configuration

**Figures (all 300 DPI PNG):**
- `Figure1.png` - Conditional means and reduced-form shocks
- `Figure2.png` - Structural shocks time series
- `Figure3.png` - Good vs bad news variance decomposition
- `Figure4.png` - Third moments (skewness) of structural shocks
- `Figure8.png` - Variance and skewness of endogenous variables
- `Figure9.png` - Covariance between GDP and headline inflation
- `Figure_GDP_Inflation_*.png` - GDP vs inflation bivariate distributions (landscape 11×8.5)
- `Figure_Unemployment_Inflation_*.png` - Unemployment vs inflation distributions (landscape 11×8.5)
- `Figure_Taylor99_and_Probabilities.png` - Federal funds rate forecast and scenario probabilities (landscape 11×8.5)
- `Figure_Macro_Quantiles.png` - Quantile time series for all macro variables (landscape 11×8.5)

---

## Technical Details

### Model Selection

**AICc (corrected Akaike Information Criterion)** is used to select among the four specifications. The correction term `(2k² + 2k)/(n_eff - k - 1)` penalizes overfitting more heavily in small samples.

**Effective sample size:** `n_eff = T/lag` accounts for overlapping observations in the lag-h forecast setting, where consecutive observations share information.

### Bootstrap Standard Errors

When `nboot > 0`:
1. Block resampling (block size = 20) preserves time series dependence
2. Each bootstrap sample is re-estimated via maximum likelihood
3. Standard errors computed from bootstrap distribution
4. Percentiles (5th, 50th, 95th) provide confidence intervals

**Computational cost:** With `nboot = 1000`, expect 2-4 hours on a modern workstation.

### Parameter Bounds

Optimization uses the following bounds:
- Scale parameters (θ_p, θ_n): [0.01, 8-10]
- Persistence (ρ_p, ρ_n): [0.01, 0.9]
- GARCH coefficients (φ_p, φ_n): [0, 8]
- Intercepts (p0, n0): [0.01, 20]

### Numerical Stability

- Initial values: pt/nt initialized at unconditional means for first `lag` observations
- Minimum values: `splus()` function enforces pt, nt ≥ 0.01
- Optimization: Grid search initialization followed by constrained optimization (fmincon)

### Likelihood Evaluation

The likelihood uses lag-period ahead forecasts: `L(t)` depends on `p(t-lag)` and `n(t-lag)`, consistent with the four-quarter-ahead forecast horizon.

---

## Advanced Usage: Filtering Mode

After initial estimation, you can use **filtering mode** for faster updates when only data changes (not model structure):

**In `estimate_bege_model.m`, set:**
```matlab
re_estimate = false;  % Use saved parameters from outbege.mat
nboot = 0;            % No bootstrap
```

**Workflow for quarterly updates:**
```matlab
% 1. Update data_pull_spfdata.csv with new quarter
% 2. Re-run the sequence
run est_structural_shocks.m
run estimate_bege_model.m    % Uses saved parameters, fast
run calc_bege_distributions.m
```

**Note:** Filtering mode requires `outbege.mat` to exist. Since this file is not included in the replication package, you must run in re-estimation mode (`re_estimate = true`) at least once.

---

## File Reference

### Main Programs (Run in Order)

1. **`est_structural_shocks.m`** - Structural shock identification
2. **`estimate_bege_model.m`** - BEGE parameter estimation
3. **`calc_bege_distributions.m`** - Forecast distribution generation

### BEGE Estimation Functions

- `BEGEspecs_inst.m` - Estimates 4 specifications, selects optimal via AICc
- `BEGEinst_llret.m` - Negative log-likelihood evaluation
- `BEGE_filter.m` - Filters pt/nt time series given parameters
- `begepdf.m` - Evaluates BEGE probability density
- `filter_shock.m` - Wrapper for filtering a single shock (used in filtering mode)
- `find_optimal_spec.m` - Extracts optimal specification from results

### Supporting Functions

- `bege_moments.m` - Calculates variance, third moment, skewness for BEGE distribution
- `bege_sim_quad.m` - Simulates quadrivariate BEGE system
- `bivconts.m` - Computes bivariate density contours via kernel smoothing
- `fmincon2.m` - Optimization wrapper (grid search + fmincon)
- `block_resample.m` - Block bootstrap resampling for time series

### Regression/VAR Functions

- `regress_boot.m` - OLS with bootstrap standard errors
- `hmoms_4factb.m` - Estimates higher moments with bootstrap covariance
- `fit_second_moms_4d.m` - GMM objective for fitting second moments

### Utility Functions

- `RecessionDummy.m` - Generates NBER recession indicator (hardcoded dates)

---

## Computational Requirements

**MATLAB:** Requires Statistics and Machine Learning Toolbox, Optimization Toolbox

**Hardware:**
- RAM: 8+ GB recommended (for 100 million simulations in `calc_bege_distributions.m`)
- CPU: Multi-core beneficial (bootstrap can use parallel computing if configured)

**Runtime estimates (typical workstation):**
- `est_structural_shocks.m`: ~30 seconds
- `estimate_bege_model.m` with `nboot = 0`: ~5-10 minutes
- `estimate_bege_model.m` with `nboot = 1000`: ~2-4 hours
- `calc_bege_distributions.m`: ~10-20 minutes (with `nsim = 1e8`)

**For testing:** Set `nsim = 1e6` in `calc_bege_distributions.m` (line 95) to reduce runtime to ~1-2 minutes.

---

## Troubleshooting

### Error: "No existing parameters found" or "Cannot load outbege.mat"

**Cause:** `re_estimate = false` but `outbege.mat` doesn't exist.

**Solution:** Set `re_estimate = true` in `estimate_bege_model.m` (line 94).

### Error: "Undefined variable rev_sup"

**Cause:** Workspace variables from `est_structural_shocks.m` not available.

**Solution:** Run `est_structural_shocks.m` first.

### Error: "Cannot open data_pull_spfdata.csv"

**Cause:** Data file not found or incorrect filename.

**Solution:** Verify `data_pull_spfdata.csv` is in the working directory. Check the filename exactly matches (including `.csv` extension).

### Out of memory during calc_bege_distributions.m

**Cause:** `nsim = 1e8` requires substantial RAM (~4-6 GB).

**Solution:** Reduce `nsim` to `1e7` or `1e6` for testing (line 95). Figures will be less smooth but adequate for verification.

---

## Citation

If you use this code, please cite:

**Engstrom, Eric (2025).** "Soft Landing or Stagflation? Measuring Macroeconomic Risks in Real Time."

---

## Code Development

This replication package was prepared and documented with the assistance of agentic artificial intelligence software and was verified for accuracy by the author.

---

## Contact

For questions or issues, please contact the author or file an issue in the project repository.