reset;

option randseed'';

#############    1. Simulierte Marktsituationen    ##############

param B := 5000;										# Anzahl Beobachtungen

param price  {1..B} := round(Uniform(1,20),1);			# Angebotspreis

param C := 5;											# Anzahl Competitors

param comp   {1..B,c in 1..C}:= round(Uniform(1,20),1);	# Preise Competitors

param rank   {i in 1..B} := 1+card({c in 1..C:comp[i,c]<price[i]});  # Rang von price 


#display price,comp;

#############    2. Simuliertes Customer Behavior    ##############

param cc1 := 0.6;
param cc2 := 0.3;
param cc3 := 0.1;

param y      {i in 1..B} :=      if rank[i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank[i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank[i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;

param test{k in 1..5} := sum{i in 1..B:rank[i]=k} y[i] / (sum{j in 1..B:rank[j]=k} 1); display test;

				
#display y; display test;

				
#############    3. Demand Learning via Logit Regression    ##############

param M default 1;											# Anzahl erklaerende Variablen

param x {i in 1..B,m in 1..M} := if m=1 then 1              # erklaerende Variablen
                           else  if m=2 then price[i] #rank[i] 

                           else  if m=3 then (if rank[i]=1 then 1)
                           else  if m=4 then (if rank[i]=2 then 1)
                           else  if m=5 then (if rank[i]=3 then 1)
						   
                           else  if m=6 then (if rank[i]>=4 then 1)
                           else  if m=7 then (price[i] - min{c in 1..C} comp[i,c]);
						   
var   beta {m in 1..M};	                                        # zu schaetzende beta-Koeffizienten

var   P {i in 1..B} = exp(sum{m in 1..M} beta[m]*x[i,m])
                  /(1+exp(sum{m in 1..M} beta[m]*x[i,m]));      # logit W'keit

subject to logMLE {m in 1..M}: sum{i in 1..B} ( y[i] - P[i] ) * x[i,m] = 0; # logit Optimalitaetsbedingungen

var LL  = sum{i in 1..B} ( y[i]*log(P[i]) + (1-y[i])*log(1-P[i]) );  # log-likelihood Funktion



#############    4. Goodness of fit    ##############

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

param AIC0 default 0;
param LL0  default 1;
param MCF0 default 0;

var MCF = 1-LL/LL0;

let M:=1; solve; let LL0 := LL; let  MCF0 := 1-LL/LL0; let AIC0 := -2*LL0 +2*1;
display beta,AIC0,LL0,MCF0;

let M:=5; solve;
display beta,AIC,LL,MCF;
                              

param P_ori {i in 1..20} := if rank[i]=1 then cc1 else if rank[i]=2 then cc2 else if rank[i]=3 then cc3 else 0;
param P_est {i in 1..20}; for {i in 1..20} let P_est[i]:=P[i];

display P_ori, P_est;


#############    5. Preisoptimierung    ##############

param N default 1;										    # Anzahl neue Beobachtungen

param pc   {i in 1..N,c in 1..C}:= round(Uniform(1,20),1);	# Preise Competitors

set A := union {i in 1..20 by 0.1} {round(i,1)};            # set of admissible prices


param rank2 {a in A,i in 1..N} := 1+card({c in 1..C:pc[i,c]<a});  # Rang von price a in situation i

param P_ori2{a in A,i in 1..N} := if rank2[a,i]=1 then cc1 else if rank2[a,i]=2 then cc2 else if rank2[a,i]=3 then cc3 else 0;
																  # originale Verkaufsw'keiten
																  
param xa {a in A,i in 1..N,m in 1..M} :=   if m=1 then 1          # erklaerende Variablen
                           else  if m=2 then a #rank2[a,i] 

                           else  if m=3 then (if rank2[a,i]=1 then 1)
                           else  if m=4 then (if rank2[a,i]=2 then 1)
                           else  if m=5 then (if rank2[a,i]=3 then 1)
						   
                           else  if m=6 then (if rank2[a,i]>=4 then 1)
                           else  if m=7 then (a - min{c in 1..C} pc[i,c]);

param bet {m in 1..M}; for {m in 1..M} let bet[m]:=beta[m];

param P_est2 {a in A,i in 1..N} := exp(sum{m in 1..M} bet[m]*xa[a,i,m])
                               /(1+exp(sum{m in 1..M} bet[m]*xa[a,i,m]));      # logit W'keit für offer price a



param Ga   {a in A}    :=             P_est2[a,N]*a;           # expected profits for situation N
param G_opt{i in 1..N} := max{a in A} P_est2[a,i]*a;           # expected profits for situation i

param a_dd{i in 1..N} := max{a in A: P_est2[a,i]*a = G_opt[i]} a;  # best price for situation i


display P_ori2, P_est2,Ga;   # teste Güte der geschätzten Verkaufsw'keiten
display pc;                  # Marktsituation N
display G_opt,a_dd;          # bester Preis + erwarteter Gewinn für situation N


#############    6. Performance Messung    ##############

param a_rand{i in 1..N} := round(Uniform(1,20),1);     # rule-based price strategies
param a_r1  {i in 1..N} := min{c in 1..C} pc[i,c];
param a_r2  {i in 1..N} := min(min{c in 1..C:a_r1[i]<pc[i,c]} pc[i,c] ,20);
param a_r3  {i in 1..N} := min{c in 1..C:a_r2[i]<pc[i,c]} pc[i,c];
param a_r12 {i in 1..N} := if Uniform(0,1)<0.5 then a_r1[i] else a_r2[i];


let N:=5000;                # "viele" realisierte zufällige Gewinne für alle Strategien

param y_dd   {i in 1..N} :=      if rank2[a_dd[i],i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank2[a_dd[i],i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank2[a_dd[i],i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;

param y_rand {i in 1..N} :=      if rank2[a_rand[i],i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank2[a_rand[i],i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank2[a_rand[i],i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;

param y_r1 {i in 1..N} :=        if rank2[a_r1[i],i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank2[a_r1[i],i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank2[a_r1[i],i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;

param y_r2 {i in 1..N} :=        if rank2[a_r2[i],i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank2[a_r2[i],i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank2[a_r2[i],i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;

param y_r12{i in 1..N} :=        if rank2[a_r12[i],i]=1 then (if Uniform(0,1)<cc1 then 1) # Ob verkauft oder nicht                       
							else if rank2[a_r12[i],i]=2 then (if Uniform(0,1)<cc2 then 1)
							else if rank2[a_r12[i],i]=3 then (if Uniform(0,1)<cc3 then 1) else 0;



param G_dd  := 1/N*sum{i in 1..N} y_dd[i]  *a_dd[i];
param G_rand:= 1/N*sum{i in 1..N} y_rand[i]*a_rand[i];
param G_r1  := 1/N*sum{i in 1..N} y_r1[i]  *a_r1[i];
param G_r2  := 1/N*sum{i in 1..N} y_r2[i]  *a_r2[i];
param G_r12 := 1/N*sum{i in 1..N} y_r12[i] *a_r12[i];

display G_dd, G_r1, G_r12, G_r2, G_rand;

#let N:=20; display a_dd;   # wie liegen die optimierten Preise so?


end;