function  val = getLambdaRVAR(gamma, cP, K0Q_cP, K1Q_cP, OmegaInv)
% restricted VAR estimation to obtain risk prices
% mu.Q, Phi.Q given
% restriction: vec( [mu, Phi] ) = R * lambda + vec( [mu^Q, phi^Q])
% estimate lambda -- Luetkepohl (5.2.6)
% All of this part comes from  Luetkepohl's book, New Introduction To
% Multiple Time Series Analysis, part 5.2. 

N = size(K0Q_cP,1) ;
T = size(cP,1)  ;

% restriction:  beta = R*lambda + r
R = zeros(N*(N+1), sum(gamma==1)) ;
R(sub2ind(size(R),(find(gamma ==1))',(1:sum(gamma==1)))) = 1 ;
r = [ K0Q_cP', reshape((K1Q_cP + eye(N)),1,[]) ]' ;  % vec( [mu^Q, phi^Q]
Z = [ones(T-1,1), cP(1:T-1,:) ]'  ;     % regressors
Y = transpose(cP(2:T, :)) ;   % dependent variable
y = reshape(Y,[],1) ; 

% beta.hat <- kronecker(solve(tcrossprod(Z))%*%Z, diag(N)) %*% y
z = y - kron(Z',eye(N)) * r ;    %the line below eq 5.2.4
inv_cov_mat = R' * kron(Z*(Z'), OmegaInv) * R  ;
cov_mat = inv(inv_cov_mat) ;
lambda_hat =  inv_cov_mat \ R' * kron(Z, OmegaInv) * z  ; % Luetkepohl (Eq 5.2.6) first equation
lambda_full = R * lambda_hat  ;   % the difference between K_p and K_Q
lam1 = lambda_full((N+1):end,:) ;

val.lambda_hat = lambda_hat ; 
val.Lam0 = lambda_full(1:N) ; 
val.Lam1=reshape(lam1, N, N);
val.cov_mat = cov_mat ;
val.inv_cov_mat = inv_cov_mat ;
val.R = R ;

end


