%
% real1.m is an example matlab program 
% to solve a simple real estate search model
%

format bank

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

clear all
diary real1.mout
diary off
delete real1.mout
diary real1.mout

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

%format short
date
time0 = clock;

% ----- set control vars ----- %

iters  = 5700;

% ----- define parameters of the model ----- %

c    = 0.07 ;
pbar = 1.90 ;
beta = (1/1.03)^(1/12) ;

beta0 = 6;
beta1 = 1;
beta2 = 1;

% ----- define states ----- %

sq   = [1;2;3] ;
bath = [1;2;3];

% - construct grid for unobserved states by quadrature - %

[ep,eppr] = qnwnorm(3,0,0.4);
[e1,epr1] = qnwnorm(3,0,0.4);
[e2,epr2] = qnwnorm(3,0,0.4);

% ----- state space: find every permuation for the states ----- %

omega=gridmake(sq,bath,ep,e1,e2);

% ----- define transition probability matrix ----- %

temp1 = eye(3,3) ;
temp2 = eye(3,3) ;

for j = 1:3
    for i = 1:3
     P(i,:,j) = prod(gridmake(temp1(:,i),temp2(:,j),eppr,epr1,epr2)') ;
 end
end

PP = [P(:,:,1) ; P(:,:,2); P(:,:,3)] ;

clear P;
P  = PP;

for j = 1:(3*3*3-1)
    P = [P;PP];
end;

% ----- find size of state space ----- %

[states,states2] = size(omega);

% ----- compute price and listing price ----- %

price    = exp(beta0 + beta1*omega(:,1) + beta2*omega(:,2) + omega(:,3)) ;
pricebar = pbar*exp(beta0 + beta1*omega(:,1) + beta2*omega(:,2) + .4) ;

% ----- compute lost commission ----- %

for j = 1:states
  if price(j,1) >= pricebar(j,1)  
   LC(j,1) = c*price(j,1);
  else 
   LC(j,1) = 0;
 end
end

% ----- value function iteration ----- %

% - initial guess of value function - %

v0  = zeros(states,1);

% - start iterations - %

crit1 = 3;

for j = 1:iters ;
 if crit1 < 0.0001, break, end;
 
% - compute second period payoffs and policy - %

  payoff1 = (1-c)*price - pricebar + omega(:,4);
  payoff2 = -LC + omega(:,5) + beta*P*v0 ;
  payoff0 = [payoff1  payoff2 ];
  
  [newv0,policy0] = max(payoff0') ;
  
% - put matrices back to correct dimension - %
 
  newv0   = newv0';
  policy0 = policy0';
  
 iter = j;

% - compute criterion - %

 crit1 = sqrt((newv0-v0)'*(newv0-v0))/(1+sqrt(v0'*v0)) ;

% - update - %

 v0 = newv0;

end;

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

[[1:states]' omega(:,1:2), price, v0, policy0],  [crit1, iter] 

% ----- save data to file ----- %

X=[[1:states]' omega(:,1:5), price, v0, policy0] ;

save real1.dat X -ascii

comptime = etime(clock, time0)
diary off