%ASSETXQ.M, a script file, not a function.
%
%   ASSETX calculates asset prices for a recursive linear equilbrium
%   with Gaussian disturbances and exponential recursive preferences
%   (see SOLVEX for the preference specification).  The parameter
%   sig governs the `risk aversion,' and is zero in the linear 
%   quadratic case.  
%   SOLVEX and ASIMUL must be run first.  Their output must be in memory.
%   The state variable x(t+1) satisfies   x(t+1) = ao x(t) + c w(t+1)
%   where ao and c are computed in SOLVEX, the asset payoff is
%   pay*x(t), discount factor is beta, and marginal utility is
%   mc*x(t).  The matrix pay and the scalar nt must be in memory.
%   The integer nt is the number of periods over which the asset prices
%   are calculated.
%   ASSETX calculates and simulates the price of a perpetual claim to
%   the stream of payouts pay*x(t).  The marginal rate of substitution,
%   the payoff, the asset price, and the gross return on the asset
%   are returned in y for a simulation of length nt.  The program
%   also computes the prices of j-period claims on sure consumption
%   for j=1,2, and 5.  The simulated values of these are returned
%   in R1, R2, and R5, respectively.


start=clock;										
disp('Asset Pricing Calculations.')
disp('      ')
  % Number of time periods to simulate is nt, which must be input
  % before using program

pay      % The asset with payoff pay*x(t) is being priced.
disp('  ')

init=x0;      %sets vector of initial conditions

init=init';

[nas,junk]=size(pay);

numer=1;        %sets numeraire as first consumption good
disp('  ')
disp('Calculating. Please Wait.')

mn=zeros(1,nc);  mn(1,numer)=1;

mn=mn*mc;


%simulate x

rand('normal');

w=rand(t1,nw);
w=u;      %Replaces w with innovations from simulation
[ys,xs]=dlsim(ao,c,[pay;mn],zeros(nas+1,nw),w,init);

%Calculate generalized variance

nbig=max(size(ao));
[rcc,ccc]=size(c);
ve=eye(ccc);
va=eye(nbig);
aomod=(va-sig*c*((ve+sig*c'*v*c)\c'*v))*ao;

[tt1,tt2]= eig(aomod);
lammax1=max(abs(diag(tt2)));
%  check on eigenvalues of aomod

if lammax1 >= 1/sqrt(beta)
disp('warning max eigenvalue of aomod exceeds 1/sqrt(beta) ')

else
end


a11s=sqrt(beta)*aomod;

%  check on eigenvalues of (ve+sig*c'*v*c)


cmod=ve+sig*c'*v*c;

eds=.0000001;   %  tolerance
% Check if cmod is positive semi-definite and symmetric
if any(eig(cmod) < -eds) | (norm(cmod'-cmod,1)/norm(cmod,1) > eds)
        disp('ccmod must be symmetric and positive semi-definite')
end
ccmod=c*inv(cmod)*c';
ar=sqrt(beta)*ao;

v1=doublej(a11s,ccmod)*beta/(1-beta);

aprice=zeros(t1,nas);
aprice3=aprice;
ii=1;
while ii <= nas;

 pay1=pay(ii,:);

 v2=doublej(a11s',pay1'*mn);
 v3=doublej2(ar',pay1',a11s',mn');


 i=1;

 aprice1=zeros(t1,1);;
 aprice2=aprice1;;
 while i <= t1;

  aprice1(i,1)=xs(i,:)*v2*xs(i,:)';
  aprice2(i,1)=xs(i,:)*v3*xs(i,:)';
  i=i+1;

 end;

 aprice(:,ii)=aprice1;

 aprice3(:,ii)=aprice2;
 ii=ii+1;

end;

clear i;


aprice=ones(t1,1)*mn*v1*pay' + aprice;
riskpr=mn*v1*pay';     % Calculates risk premium
mus1=xs*mn';

aprice=aprice./(mus1*ones(1,nas));
aprice3=aprice3./(mus1*ones(1,nas));
[ms,ns]=size(ao);
R1=beta*mn*aomod*(xs')./(mus1');
R2=(beta^2)*mn*(aomod^2)*(xs')./(mus1');
R5=(beta^5)*mn*(aomod^5)*(xs')./(mus1');
pays=xs*pay';

ret=aprice(2:t1,:)./(aprice(1:t1-1,:)-pays(1:t1-1,:));





%%%%%%%%%%%%%%%%%%%%  correction to mrs, november 25,1992

vmod=ao'*(-v + sig*v*ccmod*v)*ao;


mrs=zeros(t1-1,1);;

j=1;

while j <= t1-1
 mun=-xs(j+1,:)*v*xs(j+1,:)';
 mud=xs(j,:)*vmod*xs(j,:)';
 mun=exp(sig*mun/2)*mus1(j+1,1);
 mud=exp(sig*mud/2)*mus1(j,1);
 mrs(j,1)=mun/mud;
 j=j+1;
end;

rho=sqrt(det(cmod));
mrs=beta*rho*mrs;
pays=pays(2:t1,:);

aprice=aprice(2:t1,:);
aprice3=aprice3(2:t1,:);
timen=2:t1;

arisk=aprice-aprice3;
ya=[mrs, pays, aprice, ret, arisk];

disp('The simulated values for the marginal rate of subsitution, ...')
disp('are in the matrix ya.  The columns of ya are in the order:')
disp('[mrs, payoff, asset prices, returns, riskpr]')
disp('  ')
disp(' risk premium is not constant ')
disp(' v2 is the moment matrix in asset formula')

%%%%%%%%%% compute returns to capital


rent=(xs*md'*gamma)./(mus1*ones(1,nk));
bigk=(xs*mk')./(mus1*ones(1,nk));
bigkd=bigk*deltak;
retk=(rent(2:t1,:)+bigkd(2:t1,:))./bigk(1:t1-1,:);
yk=[mrs rent(2:t1,:) bigk(2:t1,:) retk];

disp('The simulated values for the returns to holding capital are in retk')
disp('yk=[mrs rent cap-value returns]')

t=1:(t1-1);

fin=clock;
time=fin-start;
time

%%%    std/mean of IMRS calculations

ahat=(va-2*sig*c*((ve+2*sig*c'*v*c)\c'*v))*ao;
chat=c*inv(ve +2*sig*c'*v*c)*c';

omega=ao'*(-v+2*sig*v*c*inv(ve+2*sig*c'*v*c)*c'*v)*ao;
const=mn*chat*chat'*mn';
numer=xs(1:t1-1,:)*ahat'*mn';
numer=numer.^2;
numer=numer+const*ones(t1-1,1);
detrat=det(ve +sig*c'*v*c)/sqrt(det(ve+2*sig*c'*v*c));
odiff=omega-vmod;
orat=zeros(t1-1,1);
i=1;
while i <= t1-1
 orat(i,:)=xs(i,:)*odiff*xs(i,:)';
 orat(i,:)=exp(sig*orat(i,:));
 i=i+1;
end
numer=numer.*orat;
numer=numer*detrat*beta^2;
denom=mus1(1:t1-1).^2;
mu2=numer./denom;
musq=R1(1:t1-1).^2;
sigmasq=mu2-musq';
sigma=sigmasq.^.5;
sigmurat=sigma./R1(1:t1-1)';
