% PROGRAM levelopt.m
% This MATLAB code simulates a random walk with probability p
% of going up one unit, r of remaining at the same state, and q
% of going down one unit. It assigns a cost of being in each 
% state according to a function c state^2, and assumes there is a 
% reinitialization cost K for moving the system back to state 0
% once the state hits a distance b from the origin. 
%  Thus, this is the level-crossing
% calculation in which one wants to choose the level b at which to 
% reinitialize the system so as to minimize the long-term cost of
%         (K+ A)/B
% where A is the sum of costs in the various states until level b
% is reached, and B is the time to hit level b.
%     The user is prompted to supply
% the probabilities p, and q (if p+q <1 then the walk stays in the
% current state with a positive probability), the cost K, the number n of
% simulations for each b value used, and the range of b-values
% to use (the code partitions this range into mm equal pieces). 
% At the end, output will be a vector, x, giving the times til
% the level b is reached, a vector y giving the costs along each
% sample path associated with that b value, and rat which contains the
% cost/time ratio along each path. Each of these is of length n.  
% The mean and standard deviation of each of these vectors is saved for
% each b value in the vectors mtim,stim for x, mcost, scost for y, and
% mrat,srat for rat. Finally, the long-term cost calculated as the
% ratio of the means as indicated above is calculated.
%     Graphs of these means and std are plotted as a function of the
% b values and a graph showing the relationship between the two methods
% of calculating long-term cost/time is shown.
%
!c:
p=input('p = probability of going up one unit:     ')
q=input('q = probability of going down one unit:     ')
K=input('K = reinitialization cost:   ')
c=input('c = constant in cost function:   ')
n=input('n= number of simulations for each b:'   )
blow=input('blow=lower value for range of b to reset system to 0:   ')
bhi=input('bhi=upper value for range of b to reset system to 0:   ')
% nh will count the number of moves up, nt the number of moves down
x=[ ];
y=[ ];
rat=[ ];
mrat=[ ];
srat=[ ];
mcost=[ ];
scost=[ ];
mtim=[ ];
stim=[ ];
% mm is the number of cases for b to run
mm=10;
bval(1)=blow;
for i=2:(mm+1)
   bval(i)=blow+(i-1)*(bhi-blow)/mm;
   end;
% iterate over the mm+1 b values
bval
for ii=1:(mm+1)
bb=bval(ii);
% for each b value, make n simulations
% here the current state is xn, total time is
% ttot and total cost is totcost
for j=1:n
ttot=0;
totcost=0;
xn=0;
z=[ ];
i=0;
 while (xn < bb & xn > -bb)
   i=i+1;
   a=rand(1,1);
      if a <= p
         xn=xn+1;
      elseif a >= 1-q 
         xn=xn-1;
       else
          xn=xn;
      end;
    ttot=ttot+1;
    totcost=totcost+c*xn*xn;
  z=[z xn];
 end;
% The vector x contains the times til reaching level b
% the vector y contains the total costs for each simulation 
% the vector z contains the states of the system
% the vector rat contains the sample path cost/unit time
 rr=(K+totcost)/ttot;
 x=[x ttot];
 y=[y totcost];
 rat=[rat rr];
%plot(z,'-')
%title('plot of one sample path'),pause
end;
mtim=[mtim mean(x)];
stim=[stim std(x)];
mcost=[mcost mean(y)];
scost=[scost std(y)];
mrat=[mrat mean(rat)];
srat=[srat std(rat)];
x=[ ];
y=[ ];
rat=[ ];
end;
rate=(K+mcost)./mtim;
plot(bval,mtim,'o',bval,stim,'+')
title('mean time for various b values as o, std as +'),pause
plot(bval,mcost,'o',bval,scost,'+')
title('mean cost for various b values as o, std as +'),pause
plot(bval,mrat,'o',bval,srat,'+')
title('sample path cost/time ratio as o, std as +'),pause
rate=(K+mcost)./mtim;
plot(bval,rate,'o',bval,mrat,'+')
title('ratio of means of cost/time as o, sample path costs/time as +'),pause
plot(rate,mrat,'o')
title('ratio of means of cost/time on x-axis, sample path costs/time on y-axis'),pause