%EXAMPLE 1-B  computes maximum likelihood estimates for the example in 
%             `Cattle cycles' by Rosen, et.al.   Linear trends have been
%             removed from Rosen's data.

%             Ellen R. McGrattan, 8-3-93

%File:  lqmle_1b.m 

%Data:  drosen.mat 


%Without price included:
%-----------------------

%Fixed params:

  rho1   = 0.0;
  rho2   = 0.0;
  sigma1 = 1;
  sigma2 = 1;
  x0     = [0;17.56*ones(3,1);1;0;0;0];
  grate  = 1;
  Gam    = zeros(4);


%Parameters fixed during estimation:

  delta  = 0;
  varm(4)= 1e-5;

%Current estimates:

  alpha0 = 6.393401566561099e+01;
  alpha1 = 8.895698233504461e-01;
  beta   = 9.569060069264038e-01;
  gamma0 = 3.426552943836165e-02;
  gamma1 = 3.912198387794444e-01;
  f1     = 1.024865585100780e-03;
  f3     = 1.919574544977984e-03;
  f5     = 3.118027272086529e-07;
  f7     = 5.771003859234314e-03;
  g      = 8.246179076254289e-01;
  rho3   = 7.320099995613739e-01;
  sigma3 = 1.830677216620969e+02;
  muh    = 2.085787230445125e+01;
  mum    = 7.218552156796031e+01;
  varm(1)= 8.079616958527863e-03;
  varm(2)= 1.474563399708710e-01;
  varm(3)= 7.627552695795769e+00;


%Current likelihood value:

  L      = 4.271703840841114e+02;




%With price included:
%--------------------

%Fixed params:

  rho2   = 0.0;
  sigma2 = 1;
  x0     = [0;17.56*ones(3,1);1;0;0;0];
  grate  = 1;
  Gam    = zeros(4);


%Parameters fixed during estimation:

  beta   = 0.96;
  delta  = 0;
  varm(4)= 1e-7;


%Current estimates:

  alpha0 = 1.7128e+02;  
  alpha1 = 2.3010e+00;  
  gamma0 = 2.2095e-01;  
  gamma1 = 2.1112e+00;
  f1     = 1.0247e-03; 
  f3     = 1.9193e-03; 
  f5     = 3.8389e-05; 
  f7     = 5.7707e-03; 
  g      = 8.2469e-01; 
  rho1   = 7.2214e-01; 
  rho3   = 6.7386e-01; 
  sigma1 = 1.1362e+00; 
  sigma3 = 8.2850e+00; 
  muh    = 2.8008e+01; 
  mum    = 2.4024e+01;  
  varm(1)= 8.0892e-03;
  varm(2)= 1.4760e-01; 
  varm(3)= 7.5470e+00;  
 

%Current likelihood value:

  L      = 7.3047e+02; 



%Computation of standard errors:

  param  = [alpha0;alpha1;beta;delta;gamma0;gamma1;f1;f3;f5;f7;g;rho1;rho2;
            rho3;sigma1;sigma2;sigma3;muh;mum;grate;varm(:)];
  ind    = [1,1,0,0,1,1, 1,1,1,1, 1, 1,0,1, 1,0,1, 1,1, 0, 1,1,1,0];
  ind    = [1,1,0,0,1,1, 0,0,0,0, 1, 1,0,1, 0,0,1, 0,0, 0, 1,1,1,0];

  global ZVAR
  load drosen
  format short e

  t0 = 1;
  t1 = length(ZVAR)-2;
  [L,Lt,dL,dLt]=lqnewse(t0,t1,param,ind);

  [n,m]=size(dLt);
  sum1=0;
  for t=1:m;
    sum1=sum1+dLt(:,t)*dLt(:,t)';
  end;
  se=diag(sqrt(inv(sum1)));

  disp('Results from Maximum Likelihood Estimation')
  disp('----------------------------------------')
  disp(' ')
  disp('  [Theta, Standard Errors] ')
  disp([param(ind),se])
  disp(' ')
  fprintf('  L(Theta) = %g ',L(1,1))
  disp(' ')
  disp('  Output saved in seinfo.mat')
  save seinfo

