function val = estJSZ(Y, W, N, mats, gamma)
% ## Obtain MLE
% ## Arguments:
% ##  gamma -- risk price specification: vector with length N+N^2 indicating which elements of lambda are unrestricted (1) or restricted to zero (0)
% ## Value:
% ##  pars -- list with parameter estimates

WN = W(1:N,:) ;
J = size(Y,2) ;
cP = Y * WN'  ;

if isempty(mats) == true || size(mats,2) ~= J
      disp('estML: mats needs to be provided and have length consistent with yield data') ; 
      return
end
if isempty(W) == true || size(W,2) ~= J
      disp('estML: W needs to be provided and have dimensions consistent with yield data')
      return
end

if isempty(gamma) == true
      gamma = ones(N*(N+1),1) ;
end

% (1) estimate VAR -- just to get Sigma.hat
mdl = varm(3,1) ;
lm = estimate(mdl,cP) ;
Omega_cP = lm.Covariance ;
L = chol(Omega_cP, 'lower') ;

% numerical optimization
pars_start = getStartingValues_jsz(L,Y,WN,mats) ; 
theta_start = pars2theta_jsz(pars_start)  ;
theta = getOptim(theta_start, @objJSZ,Y,WN,mats) ;
pars = theta2pars_jsz(theta, N) ;  % lamQ, Omega

% value of likelihood
res_llk = jsz_llk(Y, WN, diag(pars.lamQ-1), [], [], [], pars.Omega_cP, mats, 1, [], 1, gamma,[],[],[],[]) ;
pars.kinfQ = res_llk.K0Q_X(res_llk.K0Q_X~=0)  ;
pars.KP_0P = res_llk.K0P_cP ;
pars.KP_PP = res_llk.K1P_cP + eye(N) ;
pars.KQ_0P = res_llk.K0Q_cP  ;
pars.KQ_PP = res_llk.K1Q_cP + eye(N) ;
pars.sige2 = res_llk.sigma_e^2 ;
pars.gamma = gamma;
pars.AcP = res_llk.AcP;
pars.BcP = res_llk.BcP;
pars.rho0_cP = res_llk.rho0_cP ; 
pars.rho1_cP = res_llk.rho1_cP ;

val = pars ;

end


  