reset;

option randseed'';

param N := 1000;							# number of observations
param K default 5;							# number of explanatory variables

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

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

param y {i in 1..N} :=                     	# dependent variable y
					  if Uniform(0,1)<exp(sum{k in 1..K} b_true[k]*x[i,k])
					              /(1+exp(sum{k in 1..K} b_true[k]*x[i,k])) then 1 else 0;


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

var P {i in 1..N} = exp(sum{k in 1..K} b[k]*x[i,k])
                /(1+exp(sum{k in 1..K} b[k]*x[i,k]));  # logit probabilities

				
###########    logit optimality conditions    ###########

subject to logMLE {k in 1..K}: sum{i in 1..N} ( y[i] - P[i] ) * x[i,k] = 0;


option solver minos;

solve;

var LL  = sum{i in 1..N} ( y[i]*log(P[i]) + (1-y[i])*log(1-P[i]) );  # Log-Likelihood Function

var AIC = -2*LL +2*K;                                   # AIC

display b,b_true,LL,AIC;                                # Output


###########    McFadden Pseudo R^2   ###########

param LLK; let LLK:=LL;									# solution as param

let K:=1;												# null-model (constant only)

solve;                                                  # solve the null-model with intercept only

param LL0; let LL0:=LL;									# solution as param

param MCF = 1-LLK/LL0;

display MCF;


end;