
% farm3.m
% A quarterly version of Rosen-Murphy-Scheinkman model
%
%  formed after discussion with Jose: note that you don't need
%  big quadratic costs of c to get positive price of beef

g=.85/4;
beta=.909^(.25);
delta=.1/4  ;

gamma0=.4;     %  food factor for a calf
gamma1=.7;     %  food factor for a yearling

lambda=0; deltah=0;
thetah=0;



alpha0=500;
alpha1=.5;

pih=sqrt(1/alpha1);     %  check this  alpha1=pih^{-2}



deltak=[eye(12) zeros(12,1)];
deltak=[zeros(1,13);deltak];
deltak(1,1)=1-delta;
deltak(1,13)=g;

thetak=zeros(13,1);
thetak(1)=1;



feed=ones(12,1);

feed(1:4)=.3*ones(4,1);

feed(5:9)=.6*ones(5,1);


f=ones(28,1)*pi;

for i=1:2:27
f(i)=.001;
end
f(2)=1/(2*f(1));

for i=4:2:26
f(i)=g*feed(i/2-1)/(2*f(i-1));
end

f(28)=1/(2*f(27));


phic=zeros(15,1);
phic(1)=1;
phic(2)=f(1);
phic(15)=-f(27);

phii=zeros(15,1);
phii(1)=1;

phig=[zeros(1,14);eye(14)];

gamma=zeros(15,13);
gamma(2,1)=f(1)*(1-delta);
gamma(2,13)=g*f(1);

for i=3:12
gamma(i,i-2)=f(2*i-3);
end


rho1=0;
rho2=0;
rho3=.6^(.25);

a22=[1 0 0 0; 0 rho1 0 0; 0 0 rho2 0; 0 0 0 rho3];

c2=[0 0 0; 1 0 0; 0 1 0; 0 0 15];
c2=[0 0 0; 0 0 0; 0 0 0; 0 0 15/4];   % special case
                                    % shutting down all variances except dem.
% The conditional standard deviations of h, m, and the demand shock
% are c2(2,1), c2(3,2), and c2(4,3), respectively.


%  The first z is h, the second is m, the third is the demand shock.

uh=[50/4 1 0 0];      %  Holding costs  (divided by 4)
um=[100 0 1 0];      %  Feeding costs   (not divided by 4)


ud=[ 0 0 0 0; f(2)*uh; f(4)*uh; f(6)*uh; f(8)*uh; f(10)*uh; f(12)*uh; ...
    f(14)*uh; f(16)*uh; f(18)*uh; f(20)*uh; f(22)*uh; f(24)*uh; ...
      f(26)*uh;f(28)*um];


ub=[pih*alpha0 0 0 pih ];



solvea

% farm4.m
%  addition to farm3.m

D=zeros(2);
R=.000001*eye(2);

sy=[sk(1,:);mc];
%
%  compute steady state

nnc=15;

steadst
x0=zs;

t1=100;
nt=100;

asimul


plot(y)


A=ao;C=c;G=sy;r=4;

%function[Ar,Cr,aa,bb,cc,dd,V1]=aggreg(A,C,G,D,R,r)
%function[Ar,Cr,aa,bb,cc,dd,V1]=aggreg(A,C,G,D,R,r)
%  Does aggregation over time for a state space model.
%  The model is of the form
%
%      x(t+1) = A*x(t) + C*w(t+1)
%        y(t) = G * x(t)
%
%   However, observations (possibly error ridden) on y are available only
%   every r periods, where r is an integer greater than one.
%   The state space model that generates the observations is
%
%      x(s+1) = Ar*x(s) + Cr* er(s+1)
%        y(s) = G * x(s) + v(s)
%      v(s+1) = D * v(s) + u(s+1)
%
%    where E u(s)*u(s)'=R, E er(s)*er(s)'=Vr, and Cr = I, and where
%    s=t*r.  The program forms Ar, Cr, and Vr and uses them as inputs to
%    INNOV3.M in order to construct an innovations representation
%    for the sampled process y(s), s=1, r+1, 2r+1,... .  The innovations
%    representation is returned in the matrices aa,bb,cc,dd.  The covariance
%    matrix of the innovations is returned in V1.  Once the
%    outputs of AGGREG are in memory, VARMA2 can be used to compute
%    an arma representation for the sampled y process.
n=max(size(A));
Vr=C*C'
for i=1:r-1
Vr=Vr+(A^i)*C*C'*(A^i)';
end
Ar=A^r;
Cr=eye(n);







 [aa,bb,cc,dd,V1]=innov3(Ar,Cr,Vr,G,D,R);
										
										



%[Ar,Cr,aa,bb,cc,dd,V1]=aggreg(ao,c,sy,D,R,4);




%VARMA2.M
%   This program accepts as inputs the outputs of VARMA.M or INNOV.M.
%   This program calculates the impulse response function to each
%   innovation in the vector autoregression associated with the innovations
%   representation computed in VARMA.  It also computes impulse response
%   with respect to the orthogonalized innovations v(t) = inv(r')*u(t),
%   where r'*r = V1 = E u(t)*u(t)'.
%   The matrices aa,bb,cc,dd,V1 must be in memory.
										
%   First get impulse response for the nonorthogonalized innovations.
se=sqrt(diag(V1));
[m,n]=size(cc);
for i=1:m
i
zz=se(i)*real(dimpulse(aa,bb,cc,dd,i,40))
plot(zz),pause
end



%   Now get impulse response for the orthogonalized innovations.
rs=chol(V1);
bbb=bb*rs';
ddd=dd*rs';
for i=1:m
i
zz=real(dimpulse(aa,bbb,cc,ddd,i,40))
plot(zz),pause
end
for i=1:m
i
[n1,d1]=ss2tf(aa,bb,cc,dd,i)
end


%  farm4

%  varma2

return



%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%


nnc=5 ;
steadst;
x0=zs;
%   Simulation for x_t
sy=sk(1,:);
y1=aimpulse(ao,c,sy,1,19);
y2=aimpulse(ao,c,sy,2,19);
y3=aimpulse(ao,c,sy,3,19);

nt=150;t1=150;
asimul

%  Simulation for y(t) = x(t) + g x(t-1) + g x(t-2)

sy=-sc;
sy(1)=sy(1)+1-delta-g;
sy(2)=sy(2)+g;
sy(3)=sy(3)+g;

asimul
plot(y),title('total stock y(t)'),pause

%sy=sc;
%asimul;
%plot(y),title('consumption')



sy=[sy; sk(1,:);sc;mc];
%sy=[sk(1,:);mc];    %  price
asimul
plot(y(:,4)),title('price of beef'),pause
D=zeros(4);
R=.0001*eye(4);

sy=[sk(1,:);mc]   ;  % stock and price

nnc=5;
spectr1

return


% Now do the habit persistent part

rho_d=0;        %  shut down the serial correlation of demand

                    %  alpha1 = 1/(pih^2)
lambda=-.7*pih;     %  alpha2 =-lambda/pih  implies .7 = alpha2
deltah=.8;
thetah=1-deltah;



solvea

nnc=5 ;
steadst;
x0=zs;
%   Simulation for x_t
sy=sk(1,:);
y1=aimpulse(ao,c,sy,1,19);
y2=aimpulse(ao,c,sy,2,19);
y3=aimpulse(ao,c,sy,3,19);

nt=150;t1=150;
asimul

%  Simulation for y(t) = x(t) + g x(t-1) + g x(t-2)

sy=-sc;
sy(1)=sy(1)+1-delta-g;
sy(2)=sy(2)+g;
sy(3)=sy(3)+g;

asimul
plot(y),title('total stock y(t)'),pause

%sy=sc;
%asimul;
%plot(y),title('consumption')



sy=[sy; sk(1,:);sc;mc];

asimul
plot(y(:,2)),title('stocks'),pause
plot(y(:,4)),title('price'),pause
D=zeros(4);
R=.0001*eye(4);

%return

nnc=5;
spectr1





