% set values for parameters, initial values and data 
setvalues3
relax=.1;
relax_t=0.1;
sumL=sum(L);
sumpens=sum(pens);

expend;  % this sets K so that things are in equilibrium (working with the distribution of L and the wage levels as an input) 
         % it also sets wage levels so that on average, they are in alignment with the return to capital (the differences between regions is maintained)
         % and computes expenditures in each region, taking account of the housing multiplier

% now find the vector of wages so that supply equals demand in each region
[optw,msg]=fsolve("xssw",ones(N-1,1));
w=[1;optw];
expend

% alternatively, do this using productivity
%%%[opta,msg]=fsolve("xssa",ones(N-1,1));
%%%A=[1;opta];

% find housing stock so that prices correspond to observed prices 
Hp=house(H);
compprice=koopsom*sum(Hp)/sum(koopsom);
H=H./(compprice./Hp);
Hp=house(H);
p=price();

% now find the theta's (unobserved preferences over regions) so that the initial distibution of people is stable.

% initialize:
theta_work=zeros(N,1);
theta_pens=zeros(N,1);
[util_work,util_pens]=util(p,Hp);

% find theta_work for the distribution of workers
eps=0.0001;
theta=theta_work;
utili=util_work;
datashare=L/sum(L);
loc=locshare(utili,theta);
z=loc-datashare;
goon=1;
while (goon)
  AA=zeros(N,N);
  for jj=1:N
    theta(jj)=theta(jj)+eps;
    loc2=locshare(utili,theta);
    theta(jj)=theta(jj)-eps;
    AA(:,jj)=(loc2-loc)/(eps);
  end
  dtheta=-inv(AA)*z;
  theta=theta+relax_t*dtheta;
  theta=theta-sum(theta)/N;
  loc=locshare(utili,theta);
  z=loc-datashare;
%  [z theta];
  if (z'*z<1e-17)
    goon=0;
  end
end
theta_work=theta;
data_work=datashare;


% find theta_pens for the distribution of pensioners
theta=theta_pens;
utili=util_pens;
datashare=pens/sum(pens);
loc=locshare(utili,theta);
z=loc-datashare;
goon=1;
while (goon)
  AA=zeros(N,N);
  for jj=1:N
    theta(jj)=theta(jj)+eps;
    loc2=locshare(utili,theta);
    theta(jj)=theta(jj)-eps;
    AA(:,jj)=(loc2-loc)/(eps);
  end
  dtheta=-inv(AA)*z;
  theta=theta+relax_t*dtheta;
  theta=theta-sum(theta)/N;
  loc=locshare(utili,theta);
  z=loc-datashare;
  if (z'*z<1e-17)
    goon=0;
  end
end
theta_pens=theta;
data_pens=datashare;
%

% now print some outcomes, set some further variables
prov
printf("theta_work theta_pens\n");
[theta_work theta_pens]
init_pensshare=data_pens;
init_workshare=data_work;
init_Hp=Hp;
init_w=w;
init_p=p;
init_inhab=L+pens;
init_expend=expenditures;
init_uwork=util_work;
init_upens=util_pens;
init_L=L;
init_pens=pens;
printf("act_work act_pens kshare\n");
[data_work data_pens K/sum(K)]
w
rentincome=sum(Hp.*H)+r*sumK;
workerincome=w+(1-eta)*rentincome/sum(L)
init_workincome=workerincome;
pensincome=eta*rentincome/sum(pens)
init_pensincome=pensincome;
check=L.*workerincome+pens.*pensincome-expenditures
shar=share()
avg_own_region=sum(diag(shar))/N