% Dependent but uncorrelated: when the linear predictor fails
% Textbook: Part I, "Linear Prediction"
%
% X ~ N(0, 0.25) and Y = 0.9*X^2 + 0.1*W, W ~ N(0,1). Y is almost a
% function of X, but Cov(X, X^2) = E[X^3] = 0, so rho = 0. The best linear
% predictor is the constant E[Y], with MSE = Var[Y] = 0.11125, while the
% quadratic predictor 0.9*X^2 reaches MSE = 0.01.

clearvars; close all
rng(0)

N = 1e5;
X = 0.5*randn(1, N);
W = randn(1, N);
Y = 0.9*X.^2 + 0.1*W;

E1 = mean(X);
E2 = mean(Y);
C12 = mean(X.*Y) - E1*E2;
V1 = var(X);
rho = C12/sqrt(V1*var(Y));

Yh_lin = E2 + C12/V1*(X - E1);
Yh_quad = 0.9*X.^2;

n = 3000;
figure
plot(X(1:n), Y(1:n), '.')
hold on
x = linspace(-2, 2, 100);
plot(x, E2 + C12/V1*(x - E1), 'r', x, 0.9*x.^2, 'k--', 'LineWidth', 1.5)
hold off
xlabel('X')
ylabel('Y')
legend('data', 'linear predictor', 'quadratic predictor', 'Location', 'north')
grid on

fprintf('rho = %.4f (theory 0)\n', rho)
fprintf('MSE linear:    theory = 0.11125, simulation = %.4f\n', mean((Y - Yh_lin).^2))
fprintf('MSE quadratic: theory = 0.01,    simulation = %.4f\n', mean((Y - Yh_quad).^2))