%  HETER.M
%  Computes allocation to household i in heterogeneous agent economy.  The
%  household maximizes
%
%    -.5 E Sum Beta ^t { (si(t)-bi(t))(si(t)-bi(t)) + li(t)^2 } | I(0)
%
%  subject to
%         si(t) = lambda * hi(t-1) + pih * ci(t)
%         hi(t) = deltah * hi(t-1) + thetah * ci(t)
%    E Sum{beta^t p0(t).ci(t)}| I0 = E Sum beta^t {w0(t).li(t)
%                            + alpha0(t).di(t)}|I(0) +voki(-1)
%           bi(t)=ubi*z(t)
%           di(t)=udi*z(t)
%
%  The solution is returned in the form of matrices sci,shi,sbi,sdi, ssi,
%  sca, sha, sha, sba, sda, ssa which give the individual i and the
%  aggregate allocations for c, h, b, d, s.
%
%  The scalar tol must be set to some small nonnegative number.
%
%  SEE ALSO:  SIMULH and SOLVEA.

%STEP 1:
%Computation of the dimension of the structural parameter matrices that the
%user has specified.  The user is told the dimension of the various
%matrices in order to help the user verify that the matrices are inputted
%correctly.

	[nb,nh] = size(lambda);
	[nh,nc] = size(thetah);
	[nz,nw] = size(c2);
	[nk,ni] = size(thetak);
	[nd,nz] = size(udi);
	nx = nh+nk+nz;



%STEP 2:
%Computation of new household technology  

	a = deltah;			%Creation of the A matrix

	b = thetah;			%Creation of the B matrix

%Scaling of the law of motion by the square root of discount factor to account
%for the discounting in the optimization problem

	a = a*sqrt(beta);

	b = b*sqrt(beta);



%STEP 3:
%Creation of the matrices Q and R for the objective function of the optimal
%linear regulator problem.  The matrices Q, R, S are first created where
%the matrix S allows for the presence of cross products between the control
%and the state in the objective function.  The control is then transformed in
%the last step to remove the cross product and to put the problem in the proper
%form of the optimal linear regulator problem.


%STEP 3(i):
%Form matrices Q, R, S where S contains cross-products

	r = lambda'*lambda+tol*eye(nh);

	q = pih'*pih;
   
	s = pih'*lambda;


%STEP 3(ii):
%Transform control to remove cross-product

	if rank(q) < max(size(q));		%First check if q is nonsingular

	error('EXECUTION STOPPED.  The 1-1 block of the R matrix is singular')

	end

	qi = inv(q);

	a = a-b*qi*s;

	r = r-s'*qi*s;



%STEP 4:
%Solve the Riccati equation using two doubling algorithms.                  
%We exploit the fact that the solution to the Riccati equation can be solved
%in parts.  

	disp('Calculating, please wait')


%STEP 4(i):
%Calculate the 1-1 block of the Riccati equation.  The result is a matrix
%lambdas such that -lambdas'*h(t-1) is the feedback part of the control law.             

	[lambdas,v11] = double(a',b',r,q);

	a0 = a-b*lambdas';             

	lambdas = lambdas+qi*s;                               

	
%Step 4(ii):
%Calculate the 1-2 block of the Riccati equation.  The result is a matrix
%ubs such that ubs*z(t) is the feedforward part of the control law.

	a12 = -b*qi*pih'*ubi;

	r12 = lambda'*ubi-lambda'*pih*qi*pih'*ubi;

	r12 = r12+a0'*v11*a12;

	v12 = doublej2(a0',r12,a22',eye(nz));

	qopt = inv(q+b'*v11*b);

	v = qopt*b'*(v12*a22+v11*a12);

	ubs = v+qi*pih'*ubi;

	a0 = a0/sqrt(beta);



%Step 5:
%Compute the wealth of the individual consumer


%Step 5(i):
%Compute generalized variance

	v = doublej2(beta*ao,c,ao,c);

	v = v*beta/(1-beta);


%Step 5(ii):
%Compute the valuation of technology shock

	sdi = [zeros(nd,nh+nk),udi];

	wd = doublej2(beta*ao',md',ao',sdi');

	wd = x0'*wd*x0+sum(diag(md*v*sdi'));


%Step 5(iii):
%Compute the valuation of household input

	wg = doublej2(beta*ao',mg',ao',sg');

	wg = x0'*wg*x0+sum(diag(mg*v*sg'));


%Step 5(iv):
%Compute the valuation of consumption goods

	mhs = beta*doublej2(beta*a0',lambdas',ao',mc')*ao;

	sc1 = thetah'*mhs-mc;
	
	as = [a0 thetah*sc1; zeros(nx,nh) ao];

	x0s = [zeros(nh,1); x0];

	mcs = [zeros(nc,nh) mc];

	scs = [-lambdas sc1];

	wc1 = doublej2(beta*as',mcs',as',scs');

	cs = [zeros(nh,nw); c];

	vs = doublej2(beta*as,cs,as,cs);
	
	vs = vs*beta/(1-beta);

	wc1 = x0s'*wc1*x0s+sum(diag(mcs*vs*scs'));

	sc2 = [zeros(nc,nh+nk) ubs];

	as = [a0 thetah*sc2; zeros(nx,nh) ao];

	x0s = [h0i; x0];

	scs = [-lambdas sc2];

	wc2 = doublej2(beta*as',mcs',as',scs');

	vs = doublej2(beta*as,cs,as,cs);
	
	vs = vs*beta/(1-beta);

	wc2 = x0s'*wc2*x0s+sum(diag(mcs*vs*scs'));


%Step 5(v):
%Compute the valuation of physical capital

	wk = k0i'*(deltak'*mk+gamma'*md)*x0;



%Step 6:
%Compute the Lagrange multiplier on the budget constraint

	mu = (wk+wd-wc2)/(wc1-wg);
	


%Step 7:
%Compute the law of motion for individual consumption, household capital,
%and services

	scs = mu*sc1+sc2;

	aoi = [a0 thetah*scs; zeros(nx,nh) ao];	%ci = cs    

	sci = [-lambdas scs];

	shi = aoi(1:nh,1:nh+nx);

	ssi = [lambda zeros(nb,nx)]+pih*sci;

	sbi = [zeros(nb,2*nh+nk),ubi];

	sdi = [zeros(nd,2*nh+nk),udi];

disp('	')
disp('The matrix aoi has been calculated for the law of motion of the')
disp('augmented state vector x0s.  This matrix satisfies')
disp('	x(t+1) = aoi*x(t) + ci*w(t+1).')
disp('	')
disp('The solution to the model is given by ci(t) = sci*x(t),')
disp('hi(t) = shi*x(t), and si(t) = ssi*x(t).')
disp('	')
disp('The matrices sbi and sdi are constructed so that bi(t) = sbi*x(t),')
disp('and di(t) = sdi*x(t).')

		
%Step 8:
%Redefine solution matrices & Lagrange Multipliers obtained from SOLVEF.M
%to make them conformable with the augmented state vector x0s

	sha = [zeros(nh,nh) sh];
	ska = [zeros(nk,nh) sk];
	sia = [zeros(ni,nh) si];
	sda = [zeros(nd,nh) sd];
	sba = [zeros(nb,nh) sb];
	sca = [zeros(nc,nh) sc];
	sga = [zeros(ng,nh) sg];
	ssa = [zeros(nb,nh) ss];

	mha = [zeros(nh,nh) mh];
	mka = [zeros(nk,nh) mk];
	msa = [zeros(nb,nh) ms];
	mga = [zeros(ng,nh) mg];
	mia = [zeros(ni,nh) mi];
	mda = [zeros(nd,nh) md];
	mca = mcs;

disp('	')
disp('You are now ready to play with this economy!')
