function res = sarma(y,z,W)

%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;

x0 = zeros(2,1);

vp = eig(W);
aa = abs(vp);
W1 = W/max(aa);
W2 = -W/min(vp);

options = optimset('MaxIter',50000,'MaxFunEvals',50000,'Display','iter','TolX',1e-6);
[rho fval exitflag output] = fminsearch('LVsarma',x0,options,n,n2,y,X,W1,W2)

lambda = rho(2)
rho = rho(1)
A = eye(n2) - rho*W1;
B = inv(eye(n2) + lambda*W2);
Ay = A*y;
BAy = B*Ay;
BX = B*X;
beta = inv(BX'*BX)*BX'*BAy

res = [rho lambda]

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*W1*y;                                              % beta rho
H14 = -(1/sigma2)*X'*B'*(W2'*B' + B*W2)*resid;                                 % beta lambda
% H21 = H12';                                                               % sigma beta
H22 = (n2/(2*(sigma2^2))) - (1/(sigma2^3))*(resid'*resid);                  % sigma sigma
H23 = -(1/(sigma2^2))*y'*W1'*B'*resid;                                       % sigma rho
H24 = -(1/(sigma2^2))*resid'*W2'*B'*resid;                                    % sigma lambda
%H31 = H13';                                                                % rho beta
%H32 = H23';                                                                % rho sigma
H33 = -trace((inv(A))*W1*(inv(A))*W1) - (1/sigma2)*y'*W1'*B'*B*W1*y;            % rho rho
H34 = -(1/sigma2)*y'*W1'*B'*(W2'*B' + B*W2)*resid;                              % rho lambda
%H41 = H14';                                                                % lambda beta
%H42 = H24';                                                                % lambda sigma
%H43 = H34';                                                                % lambda rho
H44 = -trace(B*W2*B*W2) - (1/sigma2)*resid'*(W2'*B'*W2'*B'+W2'*B'*W2'*B'+W2'*B'*B*W2)*resid; % lambda lambda

H = [H11  H12  H13  H14
     H12' H22  H23  H24
     H13' H23' H33  H34
     H14' H24' H34' H44];

Infmat = inv(-H);

variance = diag(Infmat);
std = sqrt(variance)
coef = [beta
    sigma2
    rho
    lambda];
student = abs(coef./std);
pvalue = 2*(1 - norm_cdf (student))

res.nvar = nvar;
res.beta = beta;
res.rho = rho;
res.lambda = lambda;
res.lik = llike;
res.minfo = Infmat;

%%%%%%%%%%%%%%%%%%%%%%%

rhostar = rho/max(aa);
lambdastar = -lambda/min(vp);

resstar = [rhostar lambdastar]

residstar = BAy-BX*beta;
sigma2star = (1/(n2))*(residstar'*residstar)
llikestar = -(n2/2)*log(pi) - (n2/2)*log(sigma2star) + log(det(A)) - log(det(B)) - (1/(2*sigma2star))*(residstar'*residstar);

H11star = -(1/sigma2star)*X'*B'*B*X;                                                % beta beta
H12star = -(1/((sigma2star)^2))*X'*B'*residstar;                                        % beta sigma
H13star = -(1/sigma2star)*X'*B'*B*W*y;                                              % beta rho
H14star = -(1/sigma2star)*X'*B'*(W'*B' + B*W)*residstar;                                 % beta lambda
% H21star = H12star';                                                               % sigma beta
H22star = (n2/(2*(sigma2star^2))) - (1/(sigma2star^3))*(residstar'*residstar);                  % sigma sigma
H23star = -(1/(sigma2star^2))*y'*W'*B'*residstar;                                       % sigma rho
H24star = -(1/(sigma2star^2))*residstar'*W'*B'*residstar;                                    % sigma lambda
%H31star = H13star';                                                                % rho beta
%H32star = H23star';                                                                % rho sigma
H33star = -trace((inv(A))*W*(inv(A))*W) - (1/sigma2star)*y'*W'*B'*B*W*y;            % rho rho
H34star = -(1/sigma2star)*y'*W'*B'*(W'*B' + B*W)*residstar;                              % rho lambda
%H41star = H14star';                                                                % lambda beta
%H42star = H24star';                                                                % lambda sigma
%H43star = H34star';                                                                % lambda rho
H44star = -trace(B*W*B*W) - (1/sigma2star)*residstar'*(W'*B'*W'*B'+W'*B'*W'*B'+W'*B'*B*W)*residstar; % lambda lambda

Hstar = [H11star  H12star  H13star  H14star
     H12star' H22star  H23star  H24star
     H13star' H23star' H33star  H34star
     H14star' H24star' H34star' H44star];

Infmatstar = inv(-Hstar);

variancestar = diag(Infmatstar);
stdstar = sqrt(variancestar)
coefstar = [beta
    sigma2star
    rhostar
    lambdastar];
studentstar = abs(coefstar./stdstar);
pvaluestar = 2*(1 - norm_cdf (studentstar))

resstar.nvar = nvar;
resstar.beta = beta;
resstar.rhostar = rhostar;
resstar.lambdastar = lambdastar;
resstar.likstar = llikestar;
resstar.minfostar = Infmatstar;
