reset;

option randseed'';

param N := 1000;				# number of observations
param M := 0.8*N;               # number of training data
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(-2,5);		# feature data
                                # skaliere fuer Regularisierung mit stdabw in Dimension k
								
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,10);


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


minimize OLS: sum{i in 1..M} ( sum{k in 1..K} b[k]*x[i,k] - y[i] )^2;
#minimize LR1: sum{i in 1..M} ( sum{k in 1..K} b[k]*x[i,k] - y[i] )^2 + 20*sum{k in 1..K} abs(b[k]);
#minimize LR2: sum{i in 1..M} ( sum{k in 1..K} b[k]*x[i,k] - y[i] )^2 + 20*sum{k in 1..K} b[k]^2;

option solver minos;
solve;

display b_true,b;


######  goodness of fit (in sample)

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

display y_mean, variance, stddev;

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

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

var R2IS = 1-variance_new/variance;

display R2IS;


######  goodness of fit (out of sample)

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

display y_mean, y_mean2, stddev, stddev2;


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

var variance_new2 = 1/(N-M)*sum{i in M+1..N} (y[i]-fit2[i])^2;

var R2OS = 1-variance_new2/variance2;

display R2OS;