function res = sarmadouble(y,z,W,dcan,dus)

%output.algorithm The algorithm used
%output.funcCount The number of function evaluations
%output.iterationsThe number of iterations take

n2 = length(y);
n = sqrt(n2);
[m nvar] = size(z);

X = z;

Wcan = diag(dcan)*W;
Wus = diag(dus)*W;

vpcan = abs(eig(Wcan));
aacan = vpcan;
bbcan = vpcan;
for i=1:n2,
if vpcan(i)==0
    aacan(i)=-100;
    bbcan(i)=100;
end;
end;
Wcan1 = Wcan/min(bbcan);
Wcan2 = Wcan/min(bbcan);

vpus = abs(eig(Wus));
aaus = vpus;
bbus = vpus;
for i=1:n2,
if vpus(i)==0
    aaus(i)=-100;
    bbus(i)=100;
end;
end;
Wus1 = Wus/min(bbus);
Wus2 = Wus/min(bbus);

x0=zeros(4,1);

options = optimset('MaxIter',50000,'MaxFunEvals',50000,'Display','iter','TolX',1e-6);
[rho fval exitflag output] = fminsearch('LVsarmadouble',x0,options,n,n2,y,X,Wcan1,Wus1,Wcan2,Wus2)

rhocan = rho(1)/min(bbcan);
rhous = rho(2)/min(bbus);
lambdacan = rho(3)/min(bbcan);
lambdaus = rho(4)/min(bbus);

A = eye(n2) - rhocan*Wcan - rhous*Wus;
B = inv(eye(n2) + lambdacan*Wcan + lambdaus*Wus);
Ay = A*y;
BAy = B*Ay;
BX = B*X;
beta = inv(BX'*BX)*BX'*BAy;

res = [rhocan rhous lambdacan lambdaus];

resid = BAy-BX*beta;
sigma2 = (1/(n2))*(resid'*resid);

llike = -(n2/2)*log(pi) - (n2/2)*log(sigma2) + log(det(A)) - log(det(B)) - (1/(2*sigma2))*(resid'*resid);

H11 = -(1/sigma2)*X'*B'*B*X;                                                % beta beta
H12 = -(1/((sigma2)^2))*X'*B'*resid;                                        % beta sigma
H13 = -(1/sigma2)*X'*B'*B*Wcan1*y;                                           % beta rhocan
H14 = -(1/sigma2)*X'*B'*B*Wus1*y;                                            % beta rhous
H15 = -(1/sigma2)*X'*B'*(Wcan2'*B' + B*Wcan2)*resid;                           % beta lambdacan
H16 = -(1/sigma2)*X'*B'*(Wus2'*B' + B*Wus2)*resid;                             % beta lambdaus
% H21 = H12';                                                               % sigma beta
H22 = (n2/(2*(sigma2^2))) - (1/(sigma2^3))*(resid'*resid);                  % sigma sigma
H23 = -(1/(sigma2^2))*y'*Wcan1'*B'*resid;                                    % sigma rhocan
H24 = -(1/(sigma2^2))*y'*Wus1'*B'*resid;                                     % sigma rhous
H25 = -(1/(sigma2^2))*resid'*Wcan2'*B'*resid;                                 % sigma lambdacan
H26 = -(1/(sigma2^2))*resid'*Wus2'*B'*resid;                                  % sigma lambdaus
%H31 = H13';                                                                % rhocan beta
%H32 = H23';                                                                % rhocan sigma
H33 = -trace((inv(A))*Wcan1*(inv(A))*Wcan1) - (1/sigma2)*y'*Wcan1'*B'*B*Wcan1*y;% rhocan rhocan
H34 = -trace((inv(A))*Wcan1*(inv(A))*Wus1) - (1/sigma2)*y'*Wcan1'*B'*B*Wus1*y;  % rhocan rhous
H35 = -(1/sigma2)*y'*Wcan1'*B'*(Wcan2'*B' + B*Wcan2)*resid;                     % rhocan lambdacan
H36 = -(1/sigma2)*y'*Wcan1'*B'*(Wus2'*B' + B'*Wus2)*resid;                      % rhocan lambdaus
%H41 = H14';                                                                % rhous beta
%H42 = H24';                                                                % rhous sigma
%H43 = H34';                                                                % rhous rhocan
H44 = -trace((inv(A))*Wus1*(inv(A))*Wus1) - (1/sigma2)*y'*Wus1'*B'*B*Wus1*y;    % rhous rhous
H45 = -(1/sigma2)*y'*Wus1'*B'*(Wcan2'*B' + B*Wcan2)*resid;                      % rhous lambdacan
H46 = -(1/sigma2)*y'*Wus1'*B'*(Wus2'*B' + B*Wus2)*resid;                        % rhous lambdaus
%H51 = H15';                                                                % lambdacan beta
%H52 = H25';                                                                % lambdacan sigma
%H53 = H35';                                                                % lambdacan rhocan
%H54 = H45';                                                                % lambdacan rhous
H55 = -trace(B*Wcan2*B*Wcan2) - (1/sigma2)*resid'*(Wcan2'*B'*Wcan2'*B'+Wcan2'*B'*Wcan2'*B'+Wcan2'*B'*B*Wcan2)*resid; % lambdacan lambdacan
H56 = -trace(B*Wus2*B*Wcan2) - (1/sigma2)*resid'*(Wus2'*B'*Wcan2'*B'+Wcan2'*B'*Wus2'*B'+Wcan2'*B'*B*Wus2)*resid;     % lambdacan lambdaus
%H61 = H15';                                                                            % lambdaus beta
%H62 = H25';                                                                            % lambdaus sigma
%H63 = H35';                                                                            % lambdaus rhocan
%H64 = H45';                                                                            % lambdaus rhous
%H65 = H56';                                                                            % lambdaus lambdacan
H66 = -trace(B*Wus2*B*Wus2) - (1/sigma2)*resid'*(Wus2'*B'*Wus2'*B'+Wus2'*B'*Wus2'*B'+Wus2'*B'*B*Wus2)*resid; % lambdaus lambdaus


H = [H11  H12  H13  H14  H15  H16
     H12' H22  H23  H24  H25  H26
     H13' H23' H33  H34  H35  H36
     H14' H24' H34' H44  H45  H46
     H15' H25' H35' H45' H55  H56
     H16' H26' H36' H46' H56' H66];

Infmat = inv(-H);

coef = [beta
    sigma2
    rhocan
    rhous
    lambdacan
    lambdaus];
student = abs(coef./sqrt(diag(Infmat)));
pvalue = 2*(1 - norm_cdf (student))

res.nvar = nvar;
res.beta = beta;
res.pvalue = pvalue;
res.rhocan = [rhocan pvalue(123,1)];
res.rhous = [rhous pvalue(124,1)];
res.lambdacan = [lambdacan pvalue(125,1)];
res.lambdaus = [lambdaus pvalue(126,1)];
res.lik = llike;
res.minfo = Infmat;

