
% farm2.m
% A 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;
beta=.909;
delta=.1  ;

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=[(1-delta) 0 g; 1 0 0; 0 1 0];
thetak=[1 0 0 ]';


%  f5, f6 are f1,f2 in the text.  The remaining are moved up two.


f5=.001;
f6=1/(2*f5);
f1=.001;
f2=gamma0*g/(2*f1);           %  Flipped, May 16
f3=.001;
f4=gamma1*g/(2*f3);           %  Flipped, May 16
f7=.001;
f8=1/(2*f7);


phic=[1 f5 0 0 -f7]';
phii=[1 0 0 0 0]';
phig= [0 0 0 0; 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
gamma=[0 0 0; f5*(1-delta) 0 g*f5; f1 0 0; 0 f3 0; 0 0 0];

rho1=0;
rho2=0;
rho3=.6;

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];   % 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 1 0 0];      %  Holding costs
um=[100 0 1 0];      %  Feeding costs


ud=[ 0 0 0 0; f6*uh; f2*uh; f4*uh; f8*um];


ub=[pih*alpha0 0 0 pih ];

solvea

nnc=5 ;
x0 = zeros(8,1);
steadst;
x0=zs;
%   Simulation for x_t
sy=sk(1,:);
sy(2)=sy(2)+g;     %  NOTE: 2 is the position in x(t) of k_{t-1}
sy(3)=sy(3)+g;     %   Here sy is the total stock

sy=[sy;sk(1,:);sc;mc];  %  Total stock, breeding stock, consumption, price

y1=aimpulse(ao,c,sy,1,19);
y2=aimpulse(ao,c,sy,2,19);
y3=aimpulse(ao,c,sy,3,19);

nt=180;t1=180;

%  Simulation for total stock, breeding stock, consumption, price





asimul
plot([y(:,1),y(:,2)]),title('total stock and breeding stock'),pause
plot(y(:,4)),title('price of beef'),pause

yy=y;                %  Save the first simulation in yy

  return

D=zeros(4);
R=.0001*eye(4);

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

nnc=5;
 % spectr1



% 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,:);
sy(2)=sy(2)+g;     %  NOTE: 2 is the position in x(t) of k_{t-1}
sy(3)=sy(3)+g;     %   Here sy is the total stock

sy=[sy;sk(1,:);sc;mc];  %  Total stock, breeding stock, consumption, price

z1=aimpulse(ao,c,sy,1,19);
z2=aimpulse(ao,c,sy,2,19);
z3=aimpulse(ao,c,sy,3,19);

nt=80;t1=80;

%  Simulation for total stock, breeding stock, consumption, price

asimul
plot([y(:,1),y(:,2)]),title('total stock and breeding stock'),pause
plot(y(:,4)),title('price of beef'),pause

D=zeros(4);
R=.0001*eye(4);

return

nnc=5;
 %spectr1





