function irfres=impulsovecm_manu(results,T,impuls,shockspec)
%% PURPOSE: Perform IRFs to an LTV shock
%results: structure with VECM results
%T: chosen IRF horizon
%impact_unit: these are the impact responses to a one unit (100 pps) LTV shock
%impuls: a horizon's worth of LTV shocks

if results(1).coint > 0 % at least one cointegration relation

    % Main structure
    d = struct();

    %loading data

    %load('lnlr.mat'); %ln of the data in levels

    %d.data     = data;
    nlag=results(1).nlag; %Number of lags in VECM
    neqs=results(1).neqs; %Number of VECM equations
    data = results(1).data; %Data matrix
    %Setting parameters
    [nobs nvar] = size(data);
    %nend = nvar - nex; % number of endogenous variables

    %% I.- Estimating and pulling out residuals

    % i) VECM
    %  if nargin == 9
    %      x = varargin{:};
    %      results = vecm(d.data,nex,nlag,c,r,x);
    %  else
    %      results = vecm(d.data,nex,nlag,c,r);
    %  end
    %
    % We need a short run equaton for the long run exogenous variables
    % We are working with 2 data set: d2r(nex=2) and d3(nex=3)

    % ii) Short run eqs: we estimate: [nobs-(nlags+1)]

    % oil

    %      ot    = d.data(:,nvar-nex+1);          % ln of real oil price in levels
    %
    %
    %      %oil short run equation
    %
    %      dot   = tdiff(ot,1);
    %      dot   = trimr(dot,1,0);
    %      dot   = trimr(dot,nlag,0);      % diff ln of real oil price starting in the lag nlag+1
    %
    %    oilreg  = ols (dot, ones(rows(dot),1));   % as Pesaran et al, we assume real oil prices a geometric random walk
    %    oilres  = oilreg.resid;
    %


    % Pulling out parameters: norm coint space, residuals, betas dZt-j, beta Xt and alphas:

    % coint space
    cointvec = results(1).cointvec;
    x = cointvec;

    %   % 1. residuals from each equation (only the endogs)
    %  for i = 1 : nvar-nex
    %  d.err(:,i) = results(i).resid; % order is { y, q, w, fact(i),...}
    %  end
    %ADDED JAN 2019: residuals as estimated by ecm_manu1.m
    d.err = results(1).err;


    % 2. estimated coeffs of the lagged diff variables (all, endog and exog)
    for i = 1 : nvar  %I only have endogenous variables  (except ECT)
        for j = 1 : neqs*nlag %CHECK THIS LINE
            d.coeff(i,j) = results(i).beta(j)'; %reg i, coefficient on var j
            %lagcoeffs(j,i)= results(i).beta(j);
            t(j,i)= results(i).tstat(j);
        end
    end

    %disp('lagcoeffs (lagged vars x n endog):');
    % [lagcoeffs t];

    %
    %   % 3. estimated coeff of contemporanous exogenous variables
    %  for i = 1 : nend
    %      for j = nvar*nlag + 1 : nvar*nlag + nex
    %         d.coeff(i,j) = results(i).beta(j)';
    %         exogI1t(j-nvar*nlag,i)= results(i).beta(j);
    %         ext(j-nvar*nlag,i)= results(i).tstat(j);
    %     end
    %  end
    %
    %  disp('exogI1t (nex x n endog):')
    %   [exogI1t ext]

    coeff = d.coeff;
    %coeff  = [coeff; zeros(nex,rows(coeff'))]; %lagged coeffs of exogenous are zeros

    r = size(cointvec,2); %Define r: number of cointegration relations
    % 4. pulling out the coeff of adj of the long run residuals ( alphas )
    if r == 1
        for i = 1:nvar
            d.alpha(i,1) = results(i).beta(nlag*nvar+1); %THIS ONLY WORKS WITH ONE COINTEGRATION RELATION
            alpha(1,i)= results(i).beta(nlag*nvar+1);
            at(1,i)= results(i).tstat(nlag*nvar+1);
        end
    elseif r >= 2
        for i = 1:nvar
            d.alpha(i,:) = results(i).beta(nlag*nvar+1:nlag*nvar+r)';
            alpha(:,i)= results(i).beta(nlag*nvar+1:nlag*nvar+r);
            at(:,i)= results(i).tstat(nlag*nvar+1:nlag*nvar+r);
        end
    end

    % disp('alpha (num_coint_resids x n endog):');
    % [alpha at];

    %TURNED ON JAN 2019
    err    = d.err;                % residuals

    %% II.- Variance decomposition

    % Cholesky
    % rearrange err to get the selected Choleski ordering:

    %err   = [oilres(:) err(:,2) err(:,3) err(:,4) err(:,5) err(:,1)];  % :  trying diff order.... {oil, q, w, fact, y}


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

    %We decompose Omega = err(:,1 : nvar)'*err(:,1 : nvar)
    cov = (err'*err)/(rows(err)-1);
    cholcov = chol((err'*err)/(rows(err)-1))';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %   % Rearranging rows
    %     cholcovn(1,:) = cholcov(6,:); % y
    %     cholcovn(2,:) = cholcov(2,:); % q
    %     cholcovn(3,:) = cholcov(3,:); % w
    %     cholcovn(4,:) = cholcov(4,:); % oil
    %     cholcovn(5,:) = cholcov(5,:); % hw
    %     cholcovn(6,:) = cholcov(1,:); % poil
    %
    %           cholcov =  cholcovn    ;
    %
    %     % Rearranging columns
    %     cholcovn(:,1) = cholcov(:,6); % y
    %     cholcovn(:,2) = cholcov(:,2); % q
    %     cholcovn(:,3) = cholcov(:,3); % w
    %     cholcovn(:,4) = cholcov(:,4); % oil
    %     cholcovn(:,5) = cholcov(:,5); % hw
    %     cholcovn(:,6) = cholcov(:,1); % poil

    % Now rearrange Choleski matrix to get back to previous order:

    % Cholcovn is the reordered matrix

    % NOTE: Real oil price is an exog as Garrat et al. (2003)   oil a random walk

    % At this stage cholcovn is a triangular matrix with the struct shocks correlations
    % cholcovn final Choleski matrix to draw shocks from with same ordering
    % as variables in data
    % Derive impulse responses to one std
    % iterate through the effect of the transitory shock based on this:





    %%  IRF simulation
    nlag=results(1).nlag; nvar=results(1).neqs; r=results(1).coint;
    %T: Horizon over which responses tracked out (in years).time=(1:T)';


    % Declare the system variables
    DY  = zeros(T+1,nvar);   % in diffs
    Y   = zeros(T+1,nvar);   % in levels

    ECT = zeros(T+1,r);        % error correction terms

    e   = zeros(T+1,nvar);       % Shocks

    % coeff                         % Estimted param
    alpha  = d.alpha;              % Adj coeffs

    x = x(1:nvar,:); %cut off constant

    %NEXT LINE: RMAT HARD-CODED SO HAVE TO CHANGE IT BY HAND WHEN NEEDED


    if results(1).k >= 15 & results(1).k <= 17 %Cholesky robustness checks
        %Find column location of index 4 (LTV), so as to shock that index
        for i = 1 : size(results(1).chind,2)
            if results(1).chind(i) == 4
                f = i;  %Now f is the column of the LTV
            end
        end
    end

    %  Calculating IRF

    % 1. Shocks on impact
    %AGAIN, RMAT HARD-CODED IN THE LINE BELOW
    %JUNE 2022: STILL CHECK THE PARAGRAPH BELOW
    if results(1).k == 1 %3-variable VAR
        impact = cholcov(:,3)';
        impact_unit = [cholcov(:,3)./cholcov(3,3)]'; %Third column of P (for LTV shock),
        %normalized such that LTV shock = 1 (i.e. 100 pps.).
    elseif results(1).k==2 | results(1).k==3 | results(1).k==4 | results(1).k==8 | results(1).k==9 | results(1).k==10 | results(1).k==11 | results(1).k==12 | results(1).k==19 | results(1).k==20 | results(1).k==21 %4-variable VAR
        impact = cholcov(:,4)';
        impact_unit = [cholcov(:,4)./cholcov(4,4)]';
    elseif results(1).k == 5 | results(1).k == 13 | results(1).k == 14  %7-variable VAR
        impact = cholcov(:,7)';
        impact_unit = [cholcov(:,7)./cholcov(7,7)]';
    elseif results(1).k == 6 | results(1).k == 7 | results(1).k == 23 %6-variable VAR
        impact = cholcov(:,6)';
        impact_unit = [cholcov(:,6)./cholcov(6,6)]';
    elseif results(1).k == 22 %5-variable VAR
        impact = cholcov(:,5)';
        impact_unit = [cholcov(:,5)./cholcov(5,5)]';
    end

    %Alternative Cholesky ordering
    if results(1).k >= 15 & results(1).k <= 17
        impact = cholcov(:,f)'; %e.g. if LTV is the third variable, this shocks the third variable
        impact_unit = [cholcov(:,f)./cholcov(f,f)]';
    end

    if shockspec == 2 %one standard deviation shock => use first col cholcov in standard deviation terms
        e = impuls * impact; %Multiply series of impulses to structural shocks by first col cholcov
    else %Unit shocks => use first col cholcov in unit terms
        e = impuls * impact_unit;
    end


    %zer = zeros(1,nex);  % The EC term doesn't affect the exogenous

    %            d.alpha(i,1) = results(i).beta(nvar+1);
    %            alpha(1,i)= results(i).beta(nvar+1);
    DY(1,:) = e(1,:);
    Y(1,:) = DY(1,:);
    ECT(1,:) = Y(1,:) * x;

    % DY and Y
    if nlag == 1
        if r == 1
            for i= 2 : T+1
                %oilt =zeros(nvar,T+nlag);
                DY(i,:) =  [d.coeff*DY(i-1,:)']' + [d.alpha * ECT(i-1,1)]' + e(i,:);
                Y(i,:)   = Y(i-1,:) + DY(i,:);
                ECT(i,:) = Y(i,:) * x;
            end
        elseif r >= 2
            for i= 2 : T+1
                %oilt =zeros(nvar,T+nlag);
                DY(i,:) =  [d.coeff*DY(i-1,:)']' + [d.alpha * ECT(i-1,:)']' + e(i,:);
                Y(i,:)   = Y(i-1,:) + DY(i,:);
                ECT(i,:) = Y(i,:) * x;
            end
        end
    elseif nlag >= 2
        DY = [zeros(nlag-1,neqs); DY]; %Add one or more rows on top of DY for lags
        Y = [zeros(nlag-1,neqs); Y];
        ECT = [zeros(nlag-1,1); ECT];
        e = [zeros(nlag-1,neqs); e];
        for i= 2+(nlag-1) : T+1+(nlag-1)
            DYlags = nan(nlag*neqs,1);
            j = 1;
            for g = 1 : neqs
                for f = 1 : 1 : nlag
                    DYlags(j) = DY(i-f,g);
                    j = j + 1;
                end
            end

            DY(i,:) =  [d.coeff*DYlags]' + [d.alpha * ECT(i-1,1)]' + e(i,:);
            Y(i,:)   = Y(i-1,:) + DY(i,:);
            ECT(i,:) = Y(i,:) * x;
        end
    end

    %Save levels responses
    if nlag == 1
        irfres(1).respdif = DY;
        irfres(1).resplev = Y;
    elseif nlag >= 2
        irfres(1).resplev = trimr(Y,nlag-1,0); %Take first row(s) out again so shock is in row 1
    end
    %display('just a line')

    %IF CONDITION TURNED OFF JUNE 2022
    %   if shockspec <= 3 & results(1).k <= size(results(1).ltvmat,2) + 2 %4-variable with standard Cholesky
    %       % 2. plotting IRF  TURNED OFF AUGUST 2019
    % %       time=[0:T]';
    % %
    % %       figure; % We represent  100*dln(x) == %
    % %       subplot(3,2,1);plot(time,100*DY(:,1),'LineWidth',2);title('Change in LTV');ylabel('%ge point');grid on;
    % %       subplot(3,2,2);plot(time,100*DY(:,2),'LineWidth',2);title('Change in Relative house price');ylabel('%');grid on;
    % %       subplot(3,2,3);plot(time,100*DY(:,3),'LineWidth',2);title('Change in real GDP');ylabel('%');grid on;
    % %       subplot(3,2,4);plot(time,100*DY(:,4),'LineWidth',2);title('Change in Real mortgage rate');ylabel('%ge point');grid on;
    % %       subplot(3,2,5);plot(time,100*ECT(:,1),'LineWidth',2);title('EC term');ylabel('');grid on;
    % %
    % %       figure; % We represent  100*ln(x) = d in basis points
    % %       subplot(3,2,1);plot(time,100*Y(:,1),'LineWidth',2);title('LTV');ylabel('%ge point');grid on;
    % %       subplot(3,2,2);plot(time,100*Y(:,2),'LineWidth',2);title('Relative house price');ylabel('%');grid on;
    % %       subplot(3,2,3);plot(time,100*Y(:,3),'LineWidth',2);title('GDP');ylabel('%');grid on;
    % %       subplot(3,2,4);plot(time,100*Y(:,4),'LineWidth',2);title('Real mortgage rate');ylabel('%ge point');grid on;
    % %       subplot(3,2,5);plot(time,100*ECT(:,1),'LineWidth',2);title('EC term');ylabel('');grid on;
    %   end

elseif results(1).coint == 0 %No cointegration relation => a VAR


    % Main structure
    d = struct();

    %loading data

    %load('lnlr.mat'); %ln of the data in levels

    %d.data     = data;
    nlag=results(1).nlag; %Number of lags in VECM
    neqs=results(1).neqs; %Number of VECM equations
    data = results(1).data; %Data matrix
    %Setting parameters
    [nobs nvar] = size(data);
    %nend = nvar - nex; % number of endogenous variables

    %% I.- Estimating and pulling out residuals

    % i) VECM
    %  if nargin == 9
    %      x = varargin{:};
    %      results = vecm(d.data,nex,nlag,c,r,x);
    %  else
    %      results = vecm(d.data,nex,nlag,c,r);
    %  end
    %
    % We need a short run equaton for the long run exogenous variables
    % We are working with 2 data set: d2r(nex=2) and d3(nex=3)

    % ii) Short run eqs: we estimate: [nobs-(nlags+1)]

    % oil

    %      ot    = d.data(:,nvar-nex+1);          % ln of real oil price in levels
    %
    %
    %      %oil short run equation
    %
    %      dot   = tdiff(ot,1);
    %      dot   = trimr(dot,1,0);
    %      dot   = trimr(dot,nlag,0);      % diff ln of real oil price starting in the lag nlag+1
    %
    %    oilreg  = ols (dot, ones(rows(dot),1));   % as Pesaran et al, we assume real oil prices a geometric random walk
    %    oilres  = oilreg.resid;
    %


    % Pulling out parameters: 

    %   % 1. residuals from each equation (only the endogs)
    %  for i = 1 : nvar-nex
    %  d.err(:,i) = results(i).resid; % order is { y, q, w, fact(i),...}
    %  end
    %ADDED JAN 2019: residuals as estimated by ecm_manu1.m
    d.err = results(1).err;


    % 2. estimated coeffs of the lagged diff variables (all, endog and exog)
    for i = 1 : nvar  %I only have endogenous variables  (except ECT)
        for j = 1 : neqs*nlag %CHECK THIS LINE
            d.coeff(i,j) = results(i).beta(j)'; %reg i, coefficient on var j
            %lagcoeffs(j,i)= results(i).beta(j);
            t(j,i)= results(i).tstat(j);
        end
    end

    %disp('lagcoeffs (lagged vars x n endog):');
    % [lagcoeffs t];

    %
    %   % 3. estimated coeff of contemporanous exogenous variables
    %  for i = 1 : nend
    %      for j = nvar*nlag + 1 : nvar*nlag + nex
    %         d.coeff(i,j) = results(i).beta(j)';
    %         exogI1t(j-nvar*nlag,i)= results(i).beta(j);
    %         ext(j-nvar*nlag,i)= results(i).tstat(j);
    %     end
    %  end
    %
    %  disp('exogI1t (nex x n endog):')
    %   [exogI1t ext]

    coeff = d.coeff;
    %coeff  = [coeff; zeros(nex,rows(coeff'))]; %lagged coeffs of exogenous are zeros

    r = 0; %r: number of cointegration relations, which is 0 in this case



    %TURNED ON JAN 2019
    err    = d.err;                % residuals

    %% II.- Variance decomposition

    % Cholesky
    % rearrange err to get the selected Choleski ordering:

    %err   = [oilres(:) err(:,2) err(:,3) err(:,4) err(:,5) err(:,1)];  % :  trying diff order.... {oil, q, w, fact, y}


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

    %We decompose Omega = err(:,1 : nvar)'*err(:,1 : nvar)
    cov = (err'*err)/(rows(err)-1);
    cholcov = chol((err'*err)/(rows(err)-1))';
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %   % Rearranging rows
    %     cholcovn(1,:) = cholcov(6,:); % y
    %     cholcovn(2,:) = cholcov(2,:); % q
    %     cholcovn(3,:) = cholcov(3,:); % w
    %     cholcovn(4,:) = cholcov(4,:); % oil
    %     cholcovn(5,:) = cholcov(5,:); % hw
    %     cholcovn(6,:) = cholcov(1,:); % poil
    %
    %           cholcov =  cholcovn    ;
    %
    %     % Rearranging columns
    %     cholcovn(:,1) = cholcov(:,6); % y
    %     cholcovn(:,2) = cholcov(:,2); % q
    %     cholcovn(:,3) = cholcov(:,3); % w
    %     cholcovn(:,4) = cholcov(:,4); % oil
    %     cholcovn(:,5) = cholcov(:,5); % hw
    %     cholcovn(:,6) = cholcov(:,1); % poil

    % Now rearrange Choleski matrix to get back to previous order:

    % Cholcovn is the reordered matrix

    % NOTE: Real oil price is an exog as Garrat et al. (2003)   oil a random walk

    % At this stage cholcovn is a triangular matrix with the struct shocks correlations
    % cholcovn final Choleski matrix to draw shocks from with same ordering
    % as variables in data
    % Derive impulse responses to one std
    % iterate through the effect of the transitory shock based on this:





    %%  IRF simulation
    nlag=results(1).nlag; nvar=results(1).neqs; %r=results(1).coint;
    %T: Horizon over which responses tracked out (in years).time=(1:T)';


    % Declare the system variables
    DY  = zeros(T+1,nvar);   % in diffs
    Y   = zeros(T+1,nvar);   % in levels

    e   = zeros(T+1,nvar);       % Shocks

    % coeff                         % Estimted param
   

   

    %NEXT LINE: RMAT HARD-CODED SO HAVE TO CHANGE IT BY HAND WHEN NEEDED


%     if results(1).k >= 15 & results(1).k <= 17 %Cholesky robustness checks
%         %Find column location of index 4 (LTV), so as to shock that index
%         for i = 1 : size(results(1).chind,2)
%             if results(1).chind(i) == 4
%                 f = i;  %Now f is the column of the LTV
%             end
%         end
%     end

    %  Calculating IRF

    % 1. Shocks on impact
    %AGAIN, RMAT HARD-CODED IN THE LINE BELOW
    %JUNE 2022: STILL CHECK THE PARAGRAPH BELOW
%     if results(1).k == 1 %3-variable VAR
%         impact = cholcov(:,3)';
%         impact_unit = [cholcov(:,3)./cholcov(3,3)]'; %Third column of P (for LTV shock),
%         %normalized such that LTV shock = 1 (i.e. 100 pps.).
    if results(1).k == 18 
        impact = cholcov(:,4)';
        impact_unit = [cholcov(:,4)./cholcov(4,4)]';
%     elseif results(1).k == 5 | results(1).k == 13 | results(1).k == 14  %7-variable VAR
%         impact = cholcov(:,7)';
%         impact_unit = [cholcov(:,7)./cholcov(7,7)]';
%     elseif results(1).k == 6 | results(1).k == 7 %6-variable VAR
%         impact = cholcov(:,6)';
%         impact_unit = [cholcov(:,6)./cholcov(6,6)]';
    end

%     %Alternative Cholesky ordering
%     if results(1).k >= 15 & results(1).k <= 17
%         impact = cholcov(:,f)'; %e.g. if LTV is the third variable, this shocks the third variable
%         impact_unit = [cholcov(:,f)./cholcov(f,f)]';
%     end

    if shockspec == 2 %one standard deviation shock => use first col cholcov in standard deviation terms
        e = impuls * impact; %Multiply series of impulses to structural shocks by first col cholcov
    else %Unit shocks => use first col cholcov in unit terms
        e = impuls * impact_unit;
    end


    %zer = zeros(1,nex);  % The EC term doesn't affect the exogenous

    %            d.alpha(i,1) = results(i).beta(nvar+1);
    %            alpha(1,i)= results(i).beta(nvar+1);
    DY(1,:) = e(1,:);
    Y(1,:) = DY(1,:);

    % DY and Y
    if nlag == 1  %AT THIS POINT, ONLY DID VAR FOR ONE LAG
        for i= 2 : T+1
            %oilt =zeros(nvar,T+nlag);
            DY(i,:) =  [d.coeff*DY(i-1,:)']' + e(i,:);
            Y(i,:)   = Y(i-1,:) + DY(i,:);
        end
    end

    %Save levels responses
    if nlag == 1
        irfres(1).respdif = DY;
        irfres(1).resplev = Y;
    end
    %display('just a line')

    %IF CONDITION TURNED OFF JUNE 2022
    %   if shockspec <= 3 & results(1).k <= size(results(1).ltvmat,2) + 2 %4-variable with standard Cholesky
    %       % 2. plotting IRF  TURNED OFF AUGUST 2019
    % %       time=[0:T]';
    % %
    % %       figure; % We represent  100*dln(x) == %
    % %       subplot(3,2,1);plot(time,100*DY(:,1),'LineWidth',2);title('Change in LTV');ylabel('%ge point');grid on;
    % %       subplot(3,2,2);plot(time,100*DY(:,2),'LineWidth',2);title('Change in Relative house price');ylabel('%');grid on;
    % %       subplot(3,2,3);plot(time,100*DY(:,3),'LineWidth',2);title('Change in real GDP');ylabel('%');grid on;
    % %       subplot(3,2,4);plot(time,100*DY(:,4),'LineWidth',2);title('Change in Real mortgage rate');ylabel('%ge point');grid on;
    % %       subplot(3,2,5);plot(time,100*ECT(:,1),'LineWidth',2);title('EC term');ylabel('');grid on;
    % %
    % %       figure; % We represent  100*ln(x) = d in basis points
    % %       subplot(3,2,1);plot(time,100*Y(:,1),'LineWidth',2);title('LTV');ylabel('%ge point');grid on;
    % %       subplot(3,2,2);plot(time,100*Y(:,2),'LineWidth',2);title('Relative house price');ylabel('%');grid on;
    % %       subplot(3,2,3);plot(time,100*Y(:,3),'LineWidth',2);title('GDP');ylabel('%');grid on;
    % %       subplot(3,2,4);plot(time,100*Y(:,4),'LineWidth',2);title('Real mortgage rate');ylabel('%ge point');grid on;
    % %       subplot(3,2,5);plot(time,100*ECT(:,1),'LineWidth',2);title('EC term');ylabel('');grid on;
    %   end

end

return