(*Purpose of the program: simulate the distribution of wage growth in recent-increase, steady-state, and recent-decrease cases.*)

(*Written by: Yuezhou Huo*)

(*Structure of the program:
Section 1: inputs, 
  Section 2: initializations and definitions,
  Section 3: generate prices and calculate wage growth,
     Case 1: recent-increase,
     Case 2: steady-state,
     Case 3: recent-decrease,
  Section 4: plot histograms of the results.*)
(*Press Shift+Return to run the program.*)
(*All the following comments are written ABOVE the codes they describe*)

(*-----------Section 1: inputs-----------*)

(*period 0 wage*)
w0=1;
(*period 0 price*)
p0=1;
(*period 0 price level*)
pn0=1;
(*period 1 price level*)
pn=1.016;
(*lamda captures loss aversion.*)
lamda=2;
(*muA appears in the worker's utility function when the worker's surplus exceeds the firm's surplus.*)
muA=0.5;
(*muD appears in the worker's utility function when the worker's surplus is smaller than or equal to the firm's surplus.*)
muD=1;
(*sigma captures how much the worker cares about himself relative to fairness.*)
sigma=0.8;
(*cost function for the worker; e stands for effort*)
c[e_]:=e^2/4;
(*variance for period 1 real price distribution*)
eta=1;
(*number of simulations to be implemented*)
numSim=10;


(*-----------Section 2: initializations and definitions-----------*)

(*Initiate lists to hold results.*)
priceHolder={};
incWage={};
steadWage={};
decWage={};

(*Define the piecewise function used in worker's surplus function.*)
funcA[x_]=Piecewise[{{x,x≥0},{lamda*x,x<0}}];
(*Define the function for the worker's surplus.*)
funcSurplusWorker[w_,e_]:=funcA[pn*w-pn0*w0]+pn*funcA[-c[e]+c[e0]];
(*Define the function for the firm's surplus.*)
funcSurplusFirm[w_,e_]:=(pn*p*e-pn0*p0*e0)+(-pn*w+pn0*w0);
(*Define the function for the firm's profit.*)
funcFirmProf[]:=(pn*p*e)-pn*w;
(*Define the function for the worker's utility.*)
funcWorkerUtil[]:=Piecewise[{{muA*funcSurplusWorker[w]+(1-muA)*funcSurplusFirm[w]+sigma*(w0-c[e0]),funcSurplusWorker[w]>funcSurplusFirm[w]},{muD*funcSurplusWorker[w]+(1-muD)*funcSurplusFirm[w]+sigma*(w0-c[e0]),funcSurplusWorker[w]≤funcSurplusFirm[w]}}];


(*-----------Section 3: generate prices and calculate wage growth----------*)

(*This loop draws period 1 price "numSim" times and calculate wage growth for each price and for all three cases (i.e.,recent-increase, steady-state, recent-decrease)*)
For[i=1,i≤numSim,i++,
(*Draw the natural log of period 1 real price from a normal distribution with mean 0 and variance eta. Take exponential of the result as a realization of period 1 real price.*)
p=Exp[RandomVariate[NormalDistribution[0,eta]]]; 
(*Record the real price in "priceHolder"*)
AppendTo[priceHolder,p];

(*---Case 1: recent-increase---*)

(*period 0 effort.*)
e0=Extract[x/.Solve[lamda*c'[x]==p0 && x≥0,x],1] ;
(*Solve for real wage when the nominal price falls in the wage rigid region.*)
If[p0/(lamda^2)<p*pn<p0/lamda,w={w0/pn}];
(*Solve for real wage when the nominal price falls in the effort rigid region.*)
If[p0/lamda<p*pn<p0,e=e0;w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*Solve for real wage when the nominal price is above the rigid regions.*)
If[p*pn>p0,e=Extract[x/.Solve[p==lamda*c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*Solve for real wage when the nominal price is below the rigid regions.*)
If[p*pn<p0/(lamda^2),e=Extract[x/.Solve[p==c'[x]/lamda&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x],1];
(*Solve for real wage when the nominal price is exactly between the two rigid regions.*)
If[p*pn==p0/lamda,e=Extract[x/.Solve[p==c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x],1];
(*Record period 1 nominal wage in the list "incWage" only if the wage is non-negative.*)
If[w==x,w={-100}];
If[w=={-100},,w={w[[1]]}];
If[w=={-100},,AppendTo[incWage,w[[1]]*pn]];

(*The codes for steady-state and recent-decrease cases are similar to that of recent-increase case. I will only comment briefly.*)

(*---Case 2: steady-state---*)

e0=Extract[x/.Solve[c'[x]==p0 && x≥0,x],1];
(*wage rigid region*)
If[p0/lamda<p*pn<p0,w={w0/pn},e=Extract[x/.Solve[funcSurplusFirm[w,x]==funcSurplusWorker[w,x]&& x≥0,x],1]];
(*effort rigid region*)
If[p0<p*pn<lamda*p0,e=e0;w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*above rigid regions*)
If[p*pn>lamda*p0,e=Extract[x/.Solve[p==lamda*c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*below rigid regions*)
If[p*pn<p0/lamda,e=Extract[x/.Solve[p==c'[x]/lamda&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*between rigid regions*)
If[p*pn==p0,e=Extract[x/.Solve[p==c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
If[w==x,w={-100}];
If[w=={-100},,w={w[[1]]}];
If[w=={-100},,AppendTo[steadWage,w[[1]]*pn]];

(*---Case 3: recent-decrease---*)

e0=Extract[x/.Solve[c'[x]==p0*lamda&& x≥0,x],1];
(*wage rigid region*)
If[p0<p*pn<lamda*p0,w={w0/pn};e=Extract[x/.Solve[funcSurplusFirm[w,x]==funcSurplusWorker[w,x]&& x≥0,x],1]];
(*effort rigid region*)
If[lamda*p0<p*pn<lamda^2*p0,e=e0;w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*above rigid regions*)
If[p*pn>lamda^2*p0,e=Extract[x/.Solve[p==lamda*c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x]];
(*below rigid regions*)
If[p*pn<p0,e=Extract[x/.Solve[p==c'[x]/lamda&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x],1];
(*between rigid regions*)
If[p*pn==p0*lamda,e=Extract[x/.Solve[p==c'[x]&& x≥0,x],1];w=x/.Solve[funcSurplusFirm[x,e]==funcSurplusWorker[x,e]&& x≥0,x],1];
If[w==x,w={-100}];
If[w=={-100},,w={w[[1]]}];
If[w=={-100},,AppendTo[decWage,w[[1]]*pn]];
];

(*Update the wage lists from period 1 nominal wage to wage growth.*)
incWage=Log[incWage]-Log[w0];
steadWage=Log[steadWage]-Log[w0];
decWage=Log[decWage]-Log[w0];


(*-----------Section 4: plot histograms of the results-----------*)

(*Plot natural log of period 1 real prices with 50 partitions on the graph*)
Histogram[Log[priceHolder],50,AxesLabel->{log_price,frequency}]
(*Plot wage growth in recent-increase case*)
Histogram[incWage,50,AxesLabel->{wage_growth,frequency}]
(*Plot wage growth in steady-state case*)
Histogram[steadWage,50,AxesLabel->{wage_growth,frequency}]
(*Plot wage growth in recent-decrease case*)
Histogram[decWage,50,AxesLabel->{wage_growth,frequency}]

(*sample results of output:*)