function strap=bootstrap2(results,N);
%This file makes "N" bootstrapped samples using "results" from VECM estimations.
%August 2023: now also allow for possibility that it is a VAR, i.e.
%zero cointegration relations (r = 0)

%Preliminaries
nlag=results(1).nlag; %Number of lags in VECM
neqs=results(1).neqs; %Number of VECM equations
nobs=results(1).nobs; %Number of observations
r=results(1).coint; %Number of cointegration relations

strap.neqs=neqs;
strap.nobs=nobs;
strap.nlag=nlag;

%Get VECM and cointegration results

% reduced form VECM residuals
err= results(1).err; %Columns are reduced-form VECM residuals
meanerr=nan(size(err));
%Compute centered residuals
for i=1:size(err,2)
    for j=1:size(err,1)
        meanerr(j,i) = mean(err(:,i)); %Fill time series mean of an equation's error into every row
    end
end
err=err-meanerr; %now they're centered

% VECM/VAR constant term
for i=1:neqs
    con(:,i)=results(i).beta(end); %Constant from each VECM equation (stored at end of beta's), in row vector
end

if r > 0 %if there is at least one cointegration relation
    %VECM adjustment parameters to ECTs
    for i=1:neqs
        for j=neqs*nlag+1:neqs*nlag+r %In beta vectors, adjustment parameters are stored after lag coefficients
            gamma(i,j-neqs*nlag)=results(i).beta(j,:)'; %... then fetch adjustment parameters. In this case, rows indicate equations. Columns 1 through r indicate adjustment from 1st through r'th ECT.
        end
    end

    % long-run equation
    alpha=results(1).cointvec; %Cointegrating vectors in columns.
end

% coefficient matrix on lags  %WITH A VAR, THIS WORKS WITH ONE LAG - CHECK
% IF NEED MORE LAGS
for i=1:neqs
    for j=1:neqs*nlag %Lag coefficients are stored at the top of the beta's
        coeff(i,j)=results(i).beta(j,:)'; %Lag coefficients, where rows indicate equations and columns indicate lagged differenced data
        %THE j-COUNTER MAX IS NEW - CHECK
    end
end
clear i j

Xmat=nan(nobs+nlag+1,N*neqs); %Store N samples in here

%Bootstrap

%(1) First generate N sets of residuals by random sampling of residual vectors with
%replacement.


for nb=1:N %N is number of bootstrap samples
    errboot=repmat(nan,size(err)); %err are centered reduced-form VECM residuals;
    %errboot resampled residuals
    for j=1:nobs
        draw=randi(nobs); %round(0.5+nobs*rand); %One random integer between 1 and nobs
        errboot(j,:)=err(draw,:); %Draw residual vector for a particular t
    end


    % %Check distribution of draws to make sure that every residual is drawn
    % %equally often.
    % %Fine.
    % C=1e8; %number of check draws
    % checkdraw=nan(C,1);
    % for i=1:C
    %    checkdraw(i)=randi(nobs); %round(0.5+nobs*rand);
    % end
    % ndraws=zeros(nobs,1); %Vector with number of draws for each observation
    % for i=1:C
    %     for j=1:nobs
    %         if checkdraw(i)==j
    %             ndraws(j)=ndraws(j)+1;
    %         end
    %     end
    % end
    % fdraws=100*(ndraws/C); %Fraction of draws from each observation
    % figure
    % edges = [1:1:nobs];
    % histogram(checkdraw-1e-12,edges,'normalization','probability')
    % hold on
    % plot([1:nobs],repmat(mean(fdraws)/100,1,nobs),'r','Linewidth',1)

    % % %Checks to see that indeed resampled residuals.
    % %Fine.
    % check=repmat(nan,4,neqs*2);
    % figure
    % hold on
    % for i=1:neqs
    %     errser=err(:,i);
    %     errbootser=errboot(:,i);
    %     check(1,i*2-1)=min(err(:,i));  %To compare statistics of empirical and bootstrap errors
    %     check(2,i*2-1)=max(err(:,i));
    %     check(3,i*2-1)=std(err(:,i));
    %     check(4,i*2-1)=mean(err(:,i));
    %     check(1,i*2)=min(errboot(:,i));
    %     check(2,i*2)=max(errboot(:,i));
    %     check(3,i*2)=std(errboot(:,i));
    %     check(4,i*2)=mean(errboot(:,i));
    %     subplot(2,2,i)
    %     plot([1:nobs],err(:,i),[1:nobs],errboot(:,i),[1:nobs],zeros(nobs,1),'k'),legend('estimated errors','bootstrapped errors'), grid on
    %     if i==1
    %         title('VECM LTV residuals')
    %     elseif i==2
    %         title('VECM house price residuals')
    %     elseif i==3
    %         title('VECM GDP residuals')
    %     elseif i==4
    %         title('VECM mortgage rate residuals')
    %     end
    % end
    % hold off

    %(2) Then generate N sets of data
    DX=repmat(nan,nobs+nlag+1,neqs);
    X=repmat(nan,nobs+nlag+1,neqs);
    if r > 0
        ECT=repmat(nan,nobs+nlag+1,r);
    end
    e=[repmat(nan,1+nlag,neqs); errboot]; %Resampled short-run VECM residuals
    %errboot: rows are meant for 1979-2015.
    %The four above rows yield 42-row matrices/vectors.
    %Every time, the first row is meant for 1977 and the last for 2018.
    %This means: two initial values in levels. VECM with 1 lag in
    %differences => 2 lags in levels.

    if r > 0 %At least one cointegration relation
        if nlag == 1
            %Initialize at 1979 & 1980 values (with rh_ea1).
            for i=1:neqs
                X(1,i)=results(i).y(1); %Levels 1979
                X(2,i)=results(i).y(2); %Levels 1980
                DX(2,i)=X(2,i)-X(1,i); %Implied changes 1980
            end
            ECT(2,:)=[alpha'*[X(2,:)]']'; %error-correction term in 1980

            for i=3:nobs+nlag+1 %1981-2018
                DX(i,:)=con+[coeff*DX(i-1,:)']'+[gamma(:,1:r)*ECT(i-1,:)']'+e(i,:);
                X(i,:)=X(i-1,:)+DX(i,:);
                ECT(i,:)=[alpha'*[X(i,:)]']';
            end

        elseif nlag >= 2 %ALL OF THIS IS NEW - CHECK
            for i=1:neqs
                X(1,i)=results(i).y(1); %Levels 1979
                for p = 2 : nlag+1
                    X(p,i)=results(i).y(p); %Levels 1980 and up
                    DX(p,i)=X(p,i)-X(p-1,i); %Implied changes 1980 and up
                    ECT(p,:)=[alpha'*[X(p,:)]']'; %error-correction term in 1980 and up
                end
            end

            for i=nlag+2:nobs+nlag+1 %1982-2018
                DXlags = nan(nlag*neqs,1);
                j = 1;
                for g = 1 : neqs
                    for f = 1 : 1 : nlag
                        DXlags(j) = DX(i-f,g);
                        j = j + 1;
                    end
                end
                DX(i,:)=con+[coeff*DXlags]'+[gamma(:,1:r)*ECT(i-1,:)']'+e(i,:);
                X(i,:)=X(i-1,:)+DX(i,:);
                ECT(i,:)=[alpha'*[X(i,:)]']';
            end

        end

    elseif r == 0 %No cointegration relations => a VAR
        if nlag == 1 %Only do this for one lag at this stage 
            %Initialize at 1979 & 1980 values (with rh_ea1).
            for i=1:neqs
                X(1,i)=results(i).y(1); %Levels 1979
                X(2,i)=results(i).y(2); %Levels 1980
                DX(2,i)=X(2,i)-X(1,i); %Implied changes 1980
            end
   
            for i=3:nobs+nlag+1 %1981-2018
                DX(i,:)=con+[coeff*DX(i-1,:)']'+e(i,:);
                X(i,:)=X(i-1,:)+DX(i,:);        
            end
        end
    end
    %Store generated sample
    Xmat(:,(nb*neqs)-(neqs-1):nb*neqs)=X; %42 rows: from 1977 through 2018

    %     %     %See how generated data compare to actual data
    %     figure
    %     for i = 1:neqs
    %         subplot(2,2,i)
    %         plot([1:nobs+2]',[X(:,i) results(i).y]), legend('bootstrapped','actual')
    %         if i==1
    %             title('LTV')
    %         elseif i==2
    %             title('House prices')
    %         elseif i==3
    %             title('GDP')
    %         elseif i==4
    %             title('Mortgate rate')
    %         end
    %     end


end %This is end of loop capturing number of bootstrapped samples




%Xmat=trimr(Xmat,nlag,0);
strap.Xmat=Xmat; %Bootstrapped samples. 42 rows: from 1977 through 2018

%(3) Then estimate DOLS/VECM and P-T decomp/IRF/variance decomp
%for each data set: in main file

