Econometrics

% Simulated house prices (not real data): price depends on size and age, plus noise
rng(12)
n = 120;
size_m2 = 50 + 150*rand(n, 1);
age = round(60*rand(n, 1));
price = 40 + 2.5*size_m2 - 0.8*age + 25*randn(n, 1);       % thousands
T = struct('price', price, 'size_m2', size_m2, 'age', age);
mdl = fitlm(T, 'price ~ size_m2 + age')                   % coefficient table in the console
subplot(1, 2, 1)
scatter(size_m2, price, 18, age, 'filled'), colorbar
xlabel('size (m²)'), ylabel('price (k)'), title('Price by size, coloured by age')
subplot(1, 2, 2)
scatter(mdl.Fitted, mdl.Residuals, 14, 'filled'), yline(0)
xlabel('fitted'), ylabel('residual')
title(sprintf('R^2 = %.3f, +%.2f k per m² (SE %.2f)', mdl.Rsquared.Ordinary, mdl.Coefficients.Estimate(2), mdl.Coefficients.SE(2)))

% Annual mean CO2 at Mauna Loa (NOAA via github.com/datasets) — fitted with ARIMA(1,1,0) and forecast 15 years
T = readtable('https://raw.githubusercontent.com/datasets/co2-ppm/main/data/co2-annmean-mlo.csv');
y = T.Mean; yr = T.Year;
EstMdl = estimate(arima(1, 1, 0), y)
h = 15;
[yF, yMSE] = forecast(EstMdl, h, y);
yrF = yr(end) + (1:h)';
plot(yr, y, 'LineWidth', 2), hold on
plot(yrF, yF, 'r', 'LineWidth', 2)
fill([yrF; flipud(yrF)], [yF + 1.96*sqrt(yMSE); flipud(yF - 1.96*sqrt(yMSE))], 'r', 'FaceAlpha', 0.15)
hold off, legend('measured', 'forecast', '95 % interval'), xlabel('year'), ylabel('CO₂ (ppm)')
title(sprintf('ARIMA(1,1,0): %.0f ppm expected in %d', yF(end), yrF(end)))
subtitle('A statistical extrapolation, not a climate projection')

% Simulated daily returns with volatility clustering (not market data), then estimated back
rng(5)
r = simulate(garch('Constant', 2e-6, 'GARCH', 0.88, 'ARCH', 0.09), 1500);
EstMdl = estimate(garch(1, 1), r)
v = infer(EstMdl, r);                      % conditional variances
vF = forecast(EstMdl, 60, r);
subplot(2, 1, 1)
plot(100*r, 'Color', [0.6 0.6 0.6]), hold on
plot(100*1.96*sqrt(v), 'r'), plot(-100*1.96*sqrt(v), 'r'), hold off
ylabel('return (%)'), title('Returns with the model''s 95 % band')
subplot(2, 1, 2)
plot(1:1500, 100*sqrt(252*v), 1500 + (1:60), 100*sqrt(252*vF), 'LineWidth', 1.5)
ylabel('volatility (% per year)'), xlabel('day'), legend('estimated', 'forecast')
title(sprintf('Persistence α + β = %.3f', EstMdl.ARCH + EstMdl.GARCH))

% An AR(2) series: the PACF cuts off after lag 2. A random walk fails the ADF test, its changes pass.
rng(9)
y = simulate(arima('AR', [0.6 -0.3], 'Variance', 1), 400);
walk = cumsum(randn(400, 1));
subplot(2, 2, 1), plot(y), title('AR(2) series')
subplot(2, 2, 2), parcorr(y, 15)
subplot(2, 2, 3), plot(walk), title('random walk')
subplot(2, 2, 4), autocorr(walk, 15)
[~, p1] = adftest(walk); [~, p2] = adftest(diff(walk)); [~, pq] = lbqtest(y, 'Lags', 10);
sgtitle(sprintf('ADF p-value: walk %.2f, its differences %.3f · Ljung–Box on AR(2) p = %.1g', p1, p2, pq))