Methodological note – Refined copper production 2026–2031

ai-waw.com · version 1.2 · 2026-10-17 · AIWAW-2026-R01 · data contract 1.4

Contents: 1 · 2 · 3 · 4 · 5 · 6 · 7 · 8 · 9 · 10 · 11 · 12 · 13 · 14 · 15 · 16

1. Method summary and flow diagram

The forecast combines a bottom-up model – 110 mines with company guidance, 86 projects, disruptions calibrated on plan-versus-actual history, smelters and refineries, and scrap from the copper stock in use – with simple statistical models of world series. The year 2026 is set by monthly Jan–Jul data (ICSG and national agencies). Uncertainty is described by 10,000 Monte Carlo draws; interval width is calibrated on the accuracy of back-dated forecasts from 2015–2021, so it also covers errors of the model itself.

Result: world refined copper production in 2031 – median 31,806 kt (80% interval 28,732–35,324), average growth 2025–2031 +1.5% a year; mine production – 25,034 kt (+1.3% a year).

Flow diagram (to be drawn; arrows – direction of information flow):

1. Data (USGS, BGS, national agencies, companies; ICSG – current data) → historical panel 2010–2025 with class A/B/C ↓
2. Register of 110 mines and 86 projects (guidance, profiles, links, R01_*.csv assumptions) ↓
3. Monte Carlo simulation of mine output: mine paths, depletion, projects, disruptions, mines outside the register, copula ↓
4. Refining: concentrate balance, smelter capacity, country split, SX-EW, secondary refining from the stock in use ↓
5. Hierarchy nodes (countries, other countries of a region, regions, world) – coherent draws ↓
6. 2026 nowcast: prior widening, countries (EnKF), world (ICSG – importance weighting, systematic resampling) ↓
7. World-level model ensemble (equal weights) → top-down correction ↓
8. Conformal calibration: persistent level component; mass constraint (projection) ↓
9. Publication: percentiles and means from draws, CAGR, facts, scenarios, sensitivity, “100 worlds”

2. Scope and definitions

3. Data

Sources in the file: 151 (class A – official statistics and filings: 128; B – ICSG, IEA, company presentations and guidance, identities: 21; C – estimates and assumptions: 2).

3.1 Base year 2025 – source selection rules

Final national-agency data (12 months) replace the USGS “e” estimate where the definition matches (difference to USGS in the previous year ≤ 2%):

Countries outside the USGS table (“Other countries”): company data where they cover the country (|Σ assets / BGS − 1| ≤ 2%), otherwise BGS for the same year; a year without BGS – previous BGS + change in output of the same assets; BGS countries rescaled to the USGS “Other” total within ±5%.

3.2 Mine register coverage (2025, contained Cu)

3.3 Monthly 2026 data (nowcast)

3.4 Exclusions, links and class C assumptions

Register exclusions (aggregates – no double counting): 8 – M-AU-OLYMPICDAM (składnik agregatu M-AU-COPPERSA (BHP podaje wytyczne tylko dla Copper SA)); M-AU-PROMINENTHILL (składnik agregatu M-AU-COPPERSA); M-AU-CARRAPATEENA (składnik agregatu M-AU-COPPERSA); M-RS-CUKARUPEKI (składnik agregatu M-RS-ZIJINSERBIA (plan Zijin 2026 łącznie)); M-RS-BOR (składnik agregatu M-RS-ZIJINSERBIA); M-CD-KCC (od 2024 raportowana łącznie w M-CD-KCCMUTANDA); M-CD-MUTANDA (od 2024 raportowana łącznie w M-CD-KCCMUTANDA); M-AU-MTISAERNESTHENRY (agregat Glencore 2019–2020 (nieaktualny)).

Project–mine links: 32 (“replace” – profile instead of the mine path outside guidance years: 7; “add” – increment: 25).

Model constants (class C): π by stage: operating_rampup 1.00, construction 0.95, approved 0.90, feasibility 0.55, permitting 0.45, prefeasibility 0.30, pea 0.15, suspended 0.35, care_maintenance 0.30; delay multiplier by stage: operating_rampup 0.0, construction 1.0, approved 1.3, feasibility 2.0, permitting 2.2, prefeasibility 3.0, pea 3.5, suspended 2.0, care_maintenance 2.0; κ_cap = 0.95; u_max = 0.97; country rest trend ≤ 2%/yr; carry-over of the 2026 deviation to 2027: 50%; ICSG revision σ: 0.4%.

4. Mine production model

4.1 Mine path: guidance and extrapolation

For mine j and year y:

W_j(y) – company guidance midpoint (latest document, “attributable” converted to 100% where the share is known, “payable” → contained), gaps – linear interpolation; E[r] – expected performance factor from base disruption parameters (transition from plan to expected performance); g_j – depletion drift.

4.2 Depletion – empirical Bayes

b_j, s_j – log-output slope of a mature mine and its standard error (≥ 6 years); ḡ = −0.93%/yr (median of 53 mines), τ = 4.49 pp (1.4826 × MAD); SE(ḡ) = 1.2533·τ/√n = 0.77 pp – common to all mines in draw i; mines with their own slope: 38.

4.3 Projects

Ramp-up curve – hierarchical model on 14 historical ramp-ups; alternating estimation (project MAP and hyperparameters) with a Laplace correction (H_p – diagonal Hessian of the negative log-posterior at the MAP) that prevents ω from collapsing: θ̄ = (1.10, 0.46, 2.21), ω = (0.88, 0.36, 0.34), σ = 0.099; mean profile years 0–3: 15.8, 72.3, 88.9, 90.1% of capacity. Delays: n = 23 plan–start pairs, mean 0.30 years (data consistent with Poisson – r → ∞); m_stage – stage multiplier (table in 3.4). Company profiles (ramp_profile) replace the curve where given. Project plans get the full disruption factor r.

4.4 Disruptions

Simulated method of moments (12 moments of ln r: quantiles, mean, sd, tail shares) on 278 “first annual guidance – actual” pairs (one pair per mine and year – de-duplicating units reported by several owners; no smelter output); moment weights – inverse variance from a block bootstrap over units. Result: μ_ε = 0.011, σ_ε = 0.064, λ̄ = 1.216, a = 0.561, μ_x = −3.666, σ_x = 0.854, κ = 0.050 (common-component share – from the intraclass correlation 0.015); disruption allowance 3.10% vs plan, E[D] = 4.42%. Autocorrelation of ln r within units: ρ̂ = 0.120 (n = 219, t ≈ 1.78) – not significant at 5%, no AR(4.1) term. B_plan – guidance and projects; B_ext – extrapolation from actual output (already contains average disruptions, hence r/E[r]). Standard errors of the MSM parameters are not estimated (limitation – section 14).

4.5 Mines outside the register

R_c(2025) = country output − registered mines; g_c = ½ × 2015–2025 slope, bounded to [−5%, 2%]; σ_R = 10.0%/yr – median variability of Δ ln(country − registered mines with full history) 2015–2025 (21 countries; bounded to [2%, 10%]); ρ_R = 0.07 – mean correlation of country pairs within a region (21 pairs); z_reg – regional shock from the copula.

4.6 Dependence across regions

Student-t copula: R from Kendall's tau of deviations of regional log-output growth from the 5-year mean (panel of 16 years), shrunk to the mean; ρ̄ = 0.036; estimated ν = 50, bounded in the model to [4, 8] (ν = 8) as the likelihood profile on such a short sample is flat – an assumption. Sensitivity (world 2031, bottom-up model, 3,000 draws):

5. Refining model

C – Cu in concentrate (mine-model draw); φ̄ = 0.990, σ_φ = 0.022 (2010–2025), φ_max = 1.071; Cap_c – electrolytic refinery capacity (D7 register) − 2025 secondary refining, projects (new, expansions, closures), restarts and outages (year-dependent utilisation), new smelters outside the register (Manyar, Kutch – increment above output already in the base); u_c – base utilisation (3-year mean; base year for fast-growth and structural-decline countries); λ – equilibrium factor (bisection in each draw). Closing a plant with known output reduces national output by that plant's output (other plants keep their utilisation).

Groups: low – new low-cost, policy-supported smelters (China, India, Indonesia) and smelters integrated with domestic mines (domestic concentrate ≥ 50% of throughput); high – older higher-cost smelters; medium – other custom smelters. Elasticity grid on reproducing the 2023–2025 split: adopted error 10.58%, grid optimum 10.40%, version 1.0 rule 10.40% – the data do not separate the variants (differences < 2 pp), prior values adopted. Indonesian refining from concentrate ≤ φ × national concentrate (excess – to free capacity elsewhere). Country secondary shares – log-share trend (½ slope, ≤ 8%/yr) with 3% noise.

Secondary refining: lifetime distribution f(a) – lognormal mixture by sector (class C assumption):

I_t – world mine output (DS140 and panel); ρ – share of outflows reaching secondary refining, random walk from the base-year level (σ_ρ from increments of the last 15 years).

6. Model ensemble

Members: NAIVE (random walk), DRIFT (random walk with drift – mean increment ≤ 10 years), ETS (additive error, damped trend), ARIMA(1,1,0) with drift, BSTS (local trend, Gibbs sampler with FFBS) and the bottom-up model (in the backtest – the simplified BU-lite using information available at the origin). The ensemble acts on mine output and total refining (as in the backtest); the total refining correction applies to primary and secondary with the same multiplier.

Weight selection validation (leave-one-origin-out; world MAPE, h = 1–5 and h = 1):

Equal weights give the lowest out-of-sample error for mine output (3.63%) and the second lowest for refined output (3.66% vs 3.42% for the bottom-up model); out-of-sample stacking is worse (forecast combination puzzle – estimating weights on a small sample adds variance). A non-converging member is excluded from that fit (weights renormalised): 3 cases.

7. 2026 nowcast

(7.1) – ratio method within one release; σ̄ – median ratio error of countries with ≥ 4 years of history, u_k – 2025 country shares, σ_rev = 0.4%; (7.3) – prior widening, c0 = 2.62 (bottom-up coverage at h = 1); (7.4) – country EnKF (carrying 50% into 2027); (7.5) – importance weighting with the ICSG observation applied to the rest of the world (world − conditioned countries; rest σ = world σ × world x̂ / rest x̂) and systematic resampling of 10,000 from 30,000 proposals: ESS = 8,718, unique 8,155. Percentile of the observation in the prior: mine 22, refined_primary 75, refined_secondary 63.

8. Uncertainty and calibration

Base intervals in the backtest: ensemble ± z·sd of the bottom-up model (h = 1–5). A separate multiplier per stage, with out-of-sample (LOO-origin) coverage ≥ 80% and ≥ 90%; where the stage multiplier's LOO coverage falls short, the smallest safety factor f (step 0.05) that reaches it: mine output c = 2.04 (c80 = 2.04, LOO coverage 84.2%; c90 = 1.62, LOO 94.7%); refined output c = 1.94 (c80 = 1.49, LOO coverage 73.7% → safety factor f = 1.30 → 89.5%; c90 = 1.32, LOO 94.7%). For comparison, the common multiplier (version 1.1): c80 = 1.87, c90 = 1.62, bootstrap over origins 10–90% 1.46–2.05, pooled LOO coverage 76.3% / 89.5% (before calibration 55.3% and 68.4%; n = 38). c₀ (2026 prior, h = 1, 8 cases): for 90% a fallback rule – when ⌈(n+1)(1−α)⌉ > n, c = 1.25 × max E (assumption, class C). κ0, κ1 – one per draw, common to countries and stages (primary and secondary refining – world “total” scale): smooth paths, coherent draws. Multiplicative scaling on levels – zero stays zero (point mass, e.g. Panama without restart). Mass constraint (5.3): 1.17% of draws projected, violations after projection 0.

Monte Carlo error (5 seeds, full chain, 10,000 draws): standard deviation of the 2031 median across seeds – mine 41.9 kt; refined 51.9 kt.

9. Hierarchical coherence

Every draw is coherent by construction: countries and “other countries of a region” add up to regions and the world, total refining = primary + secondary (r01_publikacja.przelicz_agregaty). Adjustments (country conditioning, ensemble, widening, mass) are multipliers common to a stage's nodes or applied to bottom nodes with re-aggregation – reconciliation happens on draws (top-down), not on medians. Means (field mean) add up; percentiles do not (median of a sum ≠ sum of medians). History of other countries = region − countries shown.

10. Validation

As-of backtest (origins 2015, 2017, 2019, 2021; data vintages available at the origin); MAPE in % by horizon and mean MASE:

Validation limits: 4 independent origins (cases at different horizons of one origin are correlated); BU-lite simplifies the full model; disruption and ramp-up calibration used the full sample. Comparison with ICSG – h = 1 and common cases only (backtest table, ICSG rows with release date).

11. Scenarios and sensitivity

• Base – Median of the forecast distribution (model ensemble after the 2026 nowcast); assumptions from historical calibration. Scenarios describe alternative assumptions, not uncertainty bounds (given by the intervals). Mine disruptions as in 2010–2025 (plan-to-actual calibration). Project realisation by stage (construction 0.95; FS 0.55; PFS 0.30), negative-binomial delays. Grasberg – full capacity 2028; Kamoa-Kakula – company profile; Cobre Panamá – stockpile processing 2026–2027 (35 kt/yr), mining restart from 2028 with π 0.45; Oyu Tolgoi about 500 kt/yr from 2028. Mount Isa smelter and Townsville refinery closure in 2029; Onahama – 2027.
• Low – Tighter supply: more disruptions, weaker project delivery and longer delays. Disruption intensity ×1.4; project π ×0.7; delays ×1.5; depletion −0.5 pp/yr. Cobre Panamá not restarted; Grasberg and Kamoa-Kakula – one-year recovery slip. Smelter capacity −3%, scrap availability −5%.
• High – Smoother supply: fewer disruptions, better project delivery, faster ramp-ups. Disruption intensity ×0.7; project π ×1.25; delays ×0.7; depletion +0.4 pp/yr. Cobre Panamá – restart with π 0.9; smelter capacity +3%, scrap +5%.

Scenario = median shifted by the ratio of bottom-up medians (scenario / base, common random numbers), with the order low ≤ base ≤ high enforced. Tornado – change in the bottom-up median when one factor changes. Shapley (pick-freeze, 5 factor groups, 5 seeds):

Interpretation: a group's share of bottom-up model uncertainty (not the published one); in the file shares are rescaled to 1/c² of the published variance and the rest (1 − 1/c²: mine 76.1%, refined 73.5%) is attributed to model error (conformal calibration).

12. “100 possible worlds”

100 stratified draws: world i = the draw at quantile position (i − 0.5)/100 of the 2031 value distribution, taken from 10,000 published draws; worlds are ordered by the 2031 value and the story describes the whole world to 2031. Paths come from the same draws as the percentiles (persistent level component – no year-by-year rescaling). Drivers – actual values of the random variables in each draw; story – a template with the 2–3 factors deviating most from the median (in standard deviations), Cobre Panamá and the largest 2027–2031 disruption.

Physical path check: share of annual changes |Δ ln| > 4% in worlds ≤ 2 × share in 2010–2025 history; meta.worlds_hidden = false.

13. Comparison with ICSG

ICSG release of 2026-10-14: the “adjusted” value (after the disruption allowance and, for refining, also the concentrate availability constraint) is compared with the forecast median; TOTAL – sum of regions before adjustment. The ICSG forecast is not part of the ensemble.

Where the differences come from (each “adjusted” comparison from the latest release):

• Mine output 2026: Forecast median 23,334 kt vs ICSG “adjusted” 23,485 kt (−151 kt). 2025 base: sum of countries (USGS and national data) 23,170 kt vs 23,180 kt in ICSG data (−10 kt – definition and coverage difference). 2026 growth: forecast +0.7%, ICSG +1.3%; ICSG Jan–Jul 2026 data show −0.8% y/y, and the ICSG annual value implies +4.2% y/y in Aug–Dec.
• Mine output 2027: Forecast median 23,744 kt vs ICSG “adjusted” 24,217 kt (−473 kt). 2027 growth: forecast +1.8% y/y, ICSG +3.1%. Largest 2026→2027 changes in the forecast (draw means): Indonesia +201, China +139, DR Congo +133, Peru −92 kt; model-ensemble correction −231 kt (disruption allowance against company guidance, Grasberg and Kamoa-Kakula recovery profiles per the companies).
• Refined output 2026: Forecast median 29,246 kt vs ICSG “adjusted” 28,693 kt (+553 kt). 2025 base: sum of countries (USGS and national data) 29,089 kt vs 28,578 kt in ICSG data (+511 kt – definition and coverage difference). 2026 growth: forecast +0.5%, ICSG +0.4%; ICSG Jan–Jul 2026 data show +1.7% (primary +1.1%, secondary +4.3%) y/y, and the ICSG annual value implies −1.3% y/y in Aug–Dec.
• Refined output 2027: Forecast median 29,708 kt vs ICSG “adjusted” 29,673 kt (+35 kt). 2027 growth: forecast +1.6% y/y, ICSG +3.4%. Largest 2026→2027 changes in the forecast (draw means): China +310, Indonesia +204, Japan −186, Sweden +139 kt; model-ensemble correction −59 kt (concentrate availability, smelter capacity and secondary refining from the stock model).

Regional comparison – benchmarks table (regions not matching by area without a series link).

14. Limitations and known weaknesses

• In some countries the mine register covers less than half of output (DR Congo, China, Russia) – the country aggregate “rest” model dominates there.

• Project realisation probabilities by stage are assumptions (class C); the sample of historical ramp-ups covers only completed projects (selection bias, also in the delay distribution).

• The bottom-up backtest uses a simplified version (BU-lite); ramp-up, drift and disruption calibration used the full sample; the conformal multiplier rests on 4 independent origins.

• Smelter utilisation elasticities follow cost position (class C) – 2023–2025 data do not distinguish them.

• The country primary/secondary split in 2023–2025 is an estimate (class C); the 2026 secondary observation σ is assumed (class C).

• Chinese smelter capacity in the plant register is incomplete – a production-trend floor is used for China.

• USGS 2025 estimates (e.g. Chinese refining) may be revised; the forecast starts from that level.

• The t copula across regions is an assumption (no evidence of tail dependence); sensitivity in the self-check.

• Factors not quantified: permits and social conflicts beyond the disruption history, supply response to prices (no price data), scrap trade policy.

• Equal weights rather than stacking were chosen on the same 4 origins used to measure accuracy – mild selection (class C).

• Country rest σ: data median 12.8%/yr, capped at 10.0%/yr (variability includes revisions and coverage changes).

• The projection onto the feasible set (mass balance) trims the upper tail of primary refining – 2031 total refining p90 -9 kt; upper-tail coverage may be slightly lower.

• Standard errors of the disruption MSM parameters are not estimated; parameter uncertainty enters indirectly via conformal calibration.

• The in-use lifetime distribution (sectors) – class C assumption.

• Independent reviews A1–A5 (data, mathematics, industry, reproducibility, text) – findings and corrections made: audyt/poprawki.md, impact: audyt/poprawki-wplyw.md.

15. Reproducibility

Command: python model\run.py (options: --kalibracja, --testy, --portal, --tylko-wyniki). Environment: python 3.12.10, numpy 2.4.2, pandas 3.0.1, scipy 1.18.0, statsmodels 0.15.0, seed 20261013; exact versions – model/requirements.txt. Main seed 20261013, sensitivity seeds 20261013, 1013, 777, 4242, 9001. Calibration – out/kalibracja.json. Main run time: 818 s. Output files: out/prognoza-rafinowana.json, out/prognoza-rafinowana.md, out/metoda-rafinowana.md, out/fakty-rafinowana.json, out/backtest-rafinowana.csv, out/projekty.csv, out/kontrola_r01.json, out/mc_blad.json, out/kopula_wrazliwosc.json.

Role of AI: Claude agents (Anthropic) extracted figures from sources (with file and page, quote checks), reconciled the register, drafted code and carried out independent reviews A1–A5; all numbers are computed by deterministic code (fixed seeds).

16. Glossary, parameters, bibliography

Methodological bibliography (primary works)

• Box G.E.P., Jenkins G.M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day.
• Burgers G., van Leeuwen P.J., Evensen G. (1998). Analysis scheme in the ensemble Kalman filter. Monthly Weather Review 126, 1719–1724.
• Carter C.K., Kohn R. (1994). On Gibbs sampling for state space models. Biometrika 81(4.3), 541–553.
• Demarta S., McNeil A.J. (2005). The t copula and related copulas. International Statistical Review 73(4.1), 111–129.
• Dempster A.P., Laird N.M., Rubin D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm. JRSS B 39(4.1), 1–38.
• Efron B. (1979). Bootstrap methods: another look at the jackknife. Annals of Statistics 7(4.1), 1–26.
• Efron B., Morris C. (1975). Data analysis using Stein's estimator and its generalizations. JASA 70(350), 311–319.
• Evensen G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 99(C5), 10143–10162.
• Gardner E.S., McKenzie E. (1985). Forecasting trends in time series. Management Science 31(4.10), 1237–1246.
• Glöser S., Soulier M., Tercero Espinoza L.A. (2013). Dynamic analysis of global copper flows. Environmental Science & Technology 47(4.12), 6564–6572.
• Gordon N.J., Salmond D.J., Smith A.F.M. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F 140(4.2), 107–113.
• Hyndman R.J., Koehler A.B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting 22(4.4), 679–688.
• Hyndman R.J., Koehler A.B., Snyder R.D., Grose S. (2002). A state space framework for automatic forecasting using exponential smoothing methods. International Journal of Forecasting 18(4.3), 439–454.
• Kitagawa G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. JCGS 5(4.1), 1–25.
• Kong A., Liu J.S., Wong W.H. (1994). Sequential imputations and Bayesian missing data problems. JASA 89(425), 278–288.
• Künsch H.R. (1989). The jackknife and the bootstrap for general stationary observations. Annals of Statistics 17(4.3), 1217–1241.
• Lei J., G'Sell M., Rinaldo A., Tibshirani R.J., Wasserman L. (2018). Distribution-free predictive inference for regression. JASA 113(523), 1094–1111.
• McFadden D. (1989). A method of simulated moments for estimation of discrete response models without numerical integration. Econometrica 57(4.5), 995–1026.
• Owen A.B. (2014). Sobol' indices and Shapley value. SIAM/ASA Journal on Uncertainty Quantification 2, 245–251.
• Romano Y., Patterson E., Candès E. (2019). Conformalized quantile regression. Advances in Neural Information Processing Systems 32.
• Scott S.L., Varian H.R. (2014). Predicting the present with Bayesian structural time series. International Journal of Mathematical Modelling and Numerical Optimisation 5(1/2), 4–23.
• Shapley L.S. (1953). A value for n-person games. In: Contributions to the Theory of Games II, Annals of Mathematics Studies 28, 307–317.
• Smith J., Wallis K.F. (2009). A simple explanation of the forecast combination puzzle. Oxford Bulletin of Economics and Statistics 71(4.3), 331–355.
• Song E., Nelson B.L., Staum J. (2016). Shapley effects for global sensitivity analysis: theory and computation. SIAM/ASA JUQ 4, 1060–1083.
• Tierney L., Kadane J.B. (1986). Accurate approximations for posterior moments and marginal densities. JASA 81(393), 82–86.
• Vovk V., Gammerman A., Shafer G. (2005). Algorithmic Learning in a Random World. Springer.
• Wickramasuriya S.L., Athanasopoulos G., Hyndman R.J. (2019). Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. JASA 114(526), 804–819.
• Wolpert D.H. (1992). Stacked generalization. Neural Networks 5(4.2), 241–259.