
function [f] = shortandlong(lr,sr,masums,sigma)
%--------------------------------------------------------------------------
% Function to calculate a combination of short and long run restrictions as
% in Bjornald and Leitemo (2009), JME.  
% Code translated from RATS code of Tom Doan. 
%--------------------------------------------------------------------------
% INPUTS: lr: long-run restrictions matrix
%         sr: short-run restrictions matrix
%         masums: inv(I-B_{1}-B_{2}-...-B_{p}) where B_{i} are the VAR coefficients
%                 of lag i
%         sigma: the estimated VAR covariance matrix
%
% EXAMPLE of lr and sr:
% lr = [
% NaN NaN NaN NaN NaN;
% NaN NaN NaN NaN NaN;
% NaN NaN NaN NaN NaN;
% NaN NaN NaN NaN 0;
% NaN NaN NaN NaN NaN;];
% 
% sr = [
% NaN 0 0 0 0;
% NaN NaN 0 0 0;
% NaN NaN NaN 0 0;
% NaN NaN NaN NaN NaN;
% NaN NaN NaN NaN NaN;];
%--------------------------------------------------------------------------
% Written by Dimitris Korobilis
% University of Glasgow
% 19/04/2013
%--------------------------------------------------------------------------

iterations = 100;
cvcrit = .000001;

n = size(sigma,1);

% Compute the number of restrictions
nr = 0;
for i=1:n
    for j=1:n
        if sr(i,j)==0
            nr = nr + 1;
        end
        if lr(i,j)==0
            nr = nr + 1;
        end
    end    
end

nfree=n^2-nr;
r = zeros(nr,n^2);
fill=0;
for i=1:n
    for j=1:n
        if sr(i,j)==0
            fill = fill + 1;
            r(fill,(j-1)*n+i) = 1.0;
        end
        if lr(i,j)==0
            fill = fill + 1;
            r(fill,(j-1)*n+1:(j-1)*n+n) = masums(i,:);
        end
    end
end
rp = null(r);
%rperp = rp;

% If we want to estimate using this procedure, we need a just-identified
% model

if nr<n*(n-1)/2 || nr>n*(n-1)/2
    error('ShortAndLong Procedure Can Only Estimate Just-Identified Models')
end

% theta is the vector of free coefficients. It's initialized to the
% closest fit to the "b" vector which is the vectorized Choleski factor
% of sigma (unless the user inputs a different initial matrix.

b = chol(sigma)'; b = b(:);

theta = (rp'*rp)\rp'*b;

% Solves the system of equations sigma-b(theta)*tr(b(theta))=0 by
% Newton's method
df = zeros(nfree,nfree);
for iters = 1:iterations
    b = reshape(rp*theta,n,n);
    gap = b*b'-sigma;
    cvcrit_new = sqrt(sum(sum(gap.^2)));
    if cvcrit_new < cvcrit
        break
    end

    % Computes the matrix of derivatives of the "nfree" components of the
    % VEC'ed objective function by the components of theta.
    for i=1:nfree
        temp=b*reshape(rp(:,i),n,n)';
        b2b = temp + temp';
        df(1:nfree,i)= trilvec(b2b);
    end
    
    % Compute the Newton's method step, and adjust theta
    delta = inv(df)*trilvec(gap);
    theta = theta - delta;
end

f = b;