clear
format long e
format compact
global zxvar
load hpdata8.dat
zxvar = hpdata8;

T     = length(zxvar);
Y     = [(zxvar(:,1)+zxvar(:,7)),zxvar(:,2),zxvar(:,3), ...
        zxvar(:,5),zxvar(:,6),zxvar(:,7),zxvar(:,8),zxvar(:,9)];

mY    = Y(1,:);
d     = exp(log(Y)-hptrend(log(Y),1600));
zxvar = d.*(ones(T,1)*(mY./d(1,:)));

param=[ ...
     5.000000000000000e-01                         0
     9.990000000000001e-01                         0
     9.999000000000000e-01                         0
     1.000000000000000e-05                         0
     4.118451329759484e-01     1.000000000000000e+00
     3.264667092210679e-01     1.000000000000000e+00
     1.395903315149792e-01     1.000000000000000e+00
     2.694731027772216e-01     1.000000000000000e+00
     2.057384114638601e-01     1.000000000000000e+00
     5.395877375562486e-02     1.000000000000000e+00
     4.888062237021212e-01     1.000000000000000e+00
     5.799121269354976e+00     1.000000000000000e+00
     2.948621273270024e-02     1.000000000000000e+00
     9.867635331803416e-01     1.000000000000000e+00
     1.000000000000000e+00                         0
    -2.183640999827882e-03     1.000000000000000e+00
     7.301015565542089e-02     1.000000000000000e+00
                         0                         0
     4.534967499297793e-02     1.000000000000000e+00
     1.522472489109530e-02     1.000000000000000e+00
     1.314894444280413e+00     1.000000000000000e+00
     3.315813426318017e-01     1.000000000000000e+00
                         0                         0
     1.841619178279952e-02     1.000000000000000e+00
    -7.771953378537914e-02     1.000000000000000e+00
     1.450123127274193e-02     1.000000000000000e+00
     9.452828620741205e-01     1.000000000000000e+00
                         0                         0
    -5.171350958233063e-02     1.000000000000000e+00
    -2.114930813811186e-02     1.000000000000000e+00
                         0                         0
                         0                         0
     9.500362185514252e-01                         0
                         0                         0
                         0                         0
    -7.836341443840548e-02     1.000000000000000e+00
     1.303132271273826e-01     1.000000000000000e+00
                         0                         0
     1.378608344969311e+00     1.000000000000000e+00
     4.058982644242970e-02     1.000000000000000e+00
     2.673719360125139e-01     1.000000000000000e+00
    -7.922908794393024e-02     1.000000000000000e+00
                         0                         0
     1.253003934317140e-01     1.000000000000000e+00
     1.250693215017485e+00     1.000000000000000e+00
    -3.534582886246749e-01     1.000000000000000e+00
    -1.819689457146593e-01     1.000000000000000e+00
                         0                         0
     6.646846557918053e-03     1.000000000000000e+00
     7.387861922856144e-02     1.000000000000000e+00
    -9.202732421586281e-03     1.000000000000000e+00
     4.359448374958788e-02     1.000000000000000e+00
                         0                         0
     4.426952430739381e-02     1.000000000000000e+00
     2.138762117146142e-02     1.000000000000000e+00
                         0                         0
                         0                         0
     4.779830451464454e-02                         0
                         0                         0
                         0                         0
     1.266737253714106e-01     1.000000000000000e+00
    -2.096646710250491e-01     1.000000000000000e+00
                         0                         0
    -4.471659730136353e-01     1.000000000000000e+00
    -3.783580329233684e-02     1.000000000000000e+00
    -4.215520418740731e-01     1.000000000000000e+00
    -1.134612037193741e-01     1.000000000000000e+00
                         0                         0
    -1.125111104327259e-01     1.000000000000000e+00
    -3.423925354330764e-01     1.000000000000000e+00
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
     2.086022119950510e-03     1.000000000000000e+00
    -3.901973215481902e-04     1.000000000000000e+00
                         0                         0
     8.351910831679403e-04     1.000000000000000e+00
     1.075949026695878e-03     1.000000000000000e+00
                         0                         0
     9.868438813617021e-03     1.000000000000000e+00
                         0                         0
     9.744547774744043e-05     1.000000000000000e+00
     1.251783696597125e-04     1.000000000000000e+00
                         0                         0
                         0                         0
     5.362462611577146e-02     1.000000000000000e+00
                         0                         0
                         0                         0
                         0                         0
                         0                         0
                         0                         0
     4.359639829223389e-03     1.000000000000000e+00
     1.115332543996732e-03     1.000000000000000e+00
                         0                         0
                         0                         0
                         0                         0
                         0                         0
     8.656802353525424e-04     1.000000000000000e+00
     2.977357743939200e+02                         0
     9.231503213691488e-01     1.000000000000000e+00
     2.953371525874999e+01                         0
     1.024018769106998e+00     1.000000000000000e+00
     1.145726454895013e+01     1.000000000000000e+00
     1.106317054791681e-02     1.000000000000000e+00
     6.318658392748530e-05     1.000000000000000e+00
     2.370033719303610e-06     1.000000000000000e+00
     6.377396318999999e-01                         0
     4.009464172343737e-01     1.000000000000000e+00
     6.729778965000002e-01                         0
     4.999709733966401e-04     1.000000000000000e+00
     5.348447471937433e-01     1.000000000000000e+00
     5.010795893323747e-03     1.000000000000000e+00
     7.998670518471689e-01     1.000000000000000e+00
     7.997984291415109e-01     1.000000000000000e+00
];

ind   = param(:,2);
param = param(:,1);

[L,Lt,dL,dLt]=hpfse(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(find(ind)),se])
disp(' ')
fprintf('  L(Theta) = %g ',L)
disp(' ')
disp('  Output saved in seinfo.mat')
save seinfo

%for checking derivatives:
  
%del = diag(max(abs(param)*1e-4,1e-8));
%j=0; 
%for i=1:length(param);
%  if ind(i);
%    j=j+1;
%    [f1,f2]=hpfse(param+del(:,i),ind);
%    [m1,m2]=hpfse(param-del(:,i),ind);
%    b(i,1) = (f1-m1)/(2*del(i,i));
%    c(i,:) = (f2-m2)'/(2*del(i,i));
%  end;
%end; 
%return
