%
% rust1.m is an example matlab program 
% to compute rust's economy
%

% ----- save output to file ----- %

clear all
diary rust1.mout
diary off
delete rust1.mout
diary rust1.mout

% ----- display date and time of computation ----- %

format bank
date
time0 = clock;

% ----- set iteration controls ----- %

iters = 190800;

% ---- set state space ---- %

states = 90;

% ------ set parameters ------ %

beta    = 0.9999;
RC      = 11.7270;

theta1  = 4.8259;
theta2  = 0.577216;
theta30 = 0.3010;
theta31 = 0.6884;
theta32 = 1 - theta30 - theta31;

scale = 0.05;
alt1 = RC;

% ----- create transition probability matrix ----- %

P1 = [ theta30 theta31 theta32 ];

P = zeros(states,states);

for j = 1:(states);
 P(j,j) = theta30;
end;

for j = 1:(states-1);
 P(j,j+1) = theta31;
end;

for j = 1:(states-2);
 P(j,j+2) = theta32;
end;

% - Note: set the last probs to add to one - %

P(states,states) = 1;
P(states-1,states) = 1-theta30;

% ----- discretize milelage ----- %

x = zeros(states,1);

for j = 1:(states-1);
 x(j+1,1) = x(j,1) + 1;
end;

% ----- make the last mileage infinite ----- %

%x(states,1) = inf;
%x(states,1) = 1000000000;

% ----- value function iteration ----- %
% ----- initial matrices ----- %

v  = zeros(states,1);
v1 = zeros(states,1);
v2 = zeros(states,1);

% ----- start iterations ------ %

crit1 = 3;

for j = 1:iters ;
 if crit1 < 0.000001, break, end;

% ----- compute value function ----- %

 for i = 1:states;
  temp   = exp(v1(i,1)) + exp(v2(i,1))+eps   ;
  w(i,1) = theta2 + log(temp) ;
 end; 

% ----- compute payoffs ----- %

 for i = 1:states;
  v1(i,1) = -scale*theta1*x(i,1)+beta*P(i,:)*w ;
  v2(i,1) = -RC-scale*theta1*x(1,1)+beta*P(1,:)*w ;
  payoff(i,:)  = [v1(i,1), v2(i,1)]  ;
 end;

 [newv,policy] = max(payoff');

 newv = newv';
 policy = policy';

 iter = j;

% ----- compute criterion ----- %

 crit1 = sqrt((newv-v)'*(newv-v))/(1+sqrt(v'*v)); 

% ----- update ----- %

 v = newv;

end;

% ----- print solutions ----- %

[v, policy, x],  [crit1, iter] 

% ----- compute probabilities of choice ----- %

for i = 1:states;
 PP1(i,1) = exp(v1(i,1))/(exp(v1(i,1))+exp(v2(i,1))) ;
 PP2(i,1) = exp(v2(i,1))/(exp(v1(i,1))+exp(v2(i,1))) ;
end;

[PP1 PP2]

% ----- simulate markov chain with rule ----- %

for j = 1:10
 sim(1,j) = ceil(38*rand(1,1));
 for i = 1:1000;
  if sim(i,j) < 38 + ceil(2*(rand(1,1)-.5))
   [temp1, temp2] = markov(P,10,sim(i,j),x');
   sim(i+1,j) = temp1(1,2) + 1;
  else
   sim(i+1,j) = 1;
  end;
 end;
end;

% ----- plot data ---- %

figure(1)
plot(sim(1:100,:))
 title('(A) Mileage of Different Buses','fontsize',16,'fontname','times')
 xlabel('Month','fontsize',12,'fontname','times')    
 ylabel('Mileage','fontsize',12,'fontname','times')
axis([0 100 0 55])
pause(3)
print -depsc2 -tiff -r300 -adobecset rust1a

% ----- finish ----- %

comptime = etime(clock, time0)

diary off