reset;

option randseed'';

param N := 10000;				# number of observations
param K := 5;					# number of explanatary variables

param x {i in 1..N,k in 1..K} := if k=1 then 1 else k^2*Uniform(-1,9);		# feature data

param b_true {k in 1..K} := Uniform(-1,2);    # unknown true impact

param y {i in 1..N} := sum{k in 1..K} b_true[k]*x[i,k] + Normal(0,1);


var b {k in 1..K};	# beta coefficients (with intercept b1)


minimize OLS: sum{i in 1..N} ( sum{k in 1..K} b[k]*x[i,k] - y[i] )^2;

option solver minos;
solve;

display b_true,b;

######  goodness of fit

param y_mean  := 1/N*sum{i in 1..N}  y[i];
param variance:= 1/N*sum{i in 1..N} (y[i]-y_mean)^2;
param stddev  := if variance>0 then variance^0.5;

display y_mean, variance, stddev;


var fit{i in 1..N} = sum{k in 1..K} b[k]*x[i,k];

var variance_new = 1/N*sum{i in 1..N} (y[i]-fit[i])^2;

var R2 = 1-variance_new/variance;

display R2;