%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:  lqnew.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(3)= 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)= 7.627552695795769e+00;


%Current likelihood value:

  L      = 4.271703840841114e+02;




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

%Fixed params:

  x0     = [0;17.56*ones(3,1);1;0;0;0];
  grate  = 1;
  Gam    = zeros(4);
  f1     = 0.0001;
  f3     = 0.0001;
  f5     = 0.0001;
  f7     = 0.0001;


%Parameters fixed during estimation:

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


%Current estimates:

  alpha0 = 2.7154e+02;      % (6.6373e+01)
  alpha1 = 2.4875e+00;      % (6.4232e-01)
  gamma0 = 3.6139e-01;      % (2.3237e+00)
  gamma1 = 1.8646e+00;      % (2.4171e+00)
  g      = 8.2790e-01;      % (2.0251e-03)
  rho1   = 6.8207e-01;      % (3.7128e-02)
  rho2   = 7.3712e-01; 
  rho3   = 6.7550e-01;      % (2.4821e-02)
  sigma1 = 1.2339e+00; 
  sigma2 = 9.6492e-01; 
  sigma3 = 8.8887e+00;      % (2.2470e+00)
  muh    = 3.6922e+01; 
  mum    = 6.2970e+01; 
  varm(2)= 1.3174e-01;      % (1.3427e-02)
  varm(3)= 7.6145e+00;      % (8.9470e-01)


%Current likelihood value:

  L      = 7.1520e+02;

%Optimization of likelihood:

  Theta0 = [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, 0,0,0,0, 1, 0,1,1, 0,1,1, 0,0, 0, 1,1,0];

  Theta  = Theta0(ind);

Theta = [...
   1.4573e+02
   1.2657e+00
   6.4734e-01
   1.7706e+00
   9.3779e-01
   8.8776e-01
   6.9864e-01
   6.8173e+00
   4.0396e+00
   2.7335e-01
   4.8246e+00
];
  global ZVAR
  load drosen
  format short e
  [x,f,g]= uncmin(Theta,'lqnew');

%Computation of standard errors:

  dose = input('Do you want to compute standard errors? (1 if yes, 0 if no) ');
  if dose;
    disp(' ')
    disp('Computing standard errors .... ')
    disp(' ')
    y      = zeros(length(Theta0),1);
    ii     = find(ind);
    y(ii)  = x;
    param  = Theta0.*(~ind')+y.*ind';

    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;
    if rank(sum1)<length(sum1);
      disp('WARNING: singularity found - check if parameters identified')
      return
    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
  end;
