Error: Index in position 1 is invalid. Array indices must be positive integers or logical values?

I'm trying to use the below described code of bayesian time-varying VAR. The complete package of the code can be found here.

I have issues in estimating the TVP-VAR model with more than one lag, my data has two variables with 268 observations and the optima lag is set at 6 as indicated by information criteria.

clear; clc;
load USdata1.csv;
Y0 = USdata1(1:3,:);
shortY = USdata1(4:end,:);
[T n] = size(shortY);
Y = reshape(shortY',T*n,1);
p = 6;        % no. of lags
q = n^2*p+n;  % dim of states
Tq = T*q; Tn = T*n; q2 = q^2; nnp1 = n*(n+1);
nloop = 11000; burnin = 1000;
    %% initialize for storage
store_Omega11 = zeros(nloop-burnin,n,n);
store_Omega22 = zeros(nloop-burnin,q);
store_beta = zeros(nloop-burnin,Tq);
    %% initialize the Markov chain
Omega11 = cov(USdata1);
Omega22 = .01*ones(q,1); % store only the diag elements
invOmega11 = Omega11\speye(n);
invOmega22 = 1./Omega22;
beta = zeros(Tq,1);
    %% prior
nu01 = n+3; S01 = eye(n);
nu02 = 3;  S02 = .005*ones(q,1);
Vbeta = 5*ones(q,1);
    %% construct and compute a few thigns
X = zeros(T,n*p);
for i=1:p
    X(:,(i-1)*n+1:i*n) = [Y0(3-i+1:end,:); shortY(1:T-i,:)];
end
bigG = SURform([ones(n*T,1) kron(X,ones(n,1))]);
H = speye(Tq,Tq) - sparse(q+1:Tq,1:(T-1)*q,ones(1,(T-1)*q),Tq,Tq);
newnu1 = nu01 + T;
newnu2 = nu02 + (T-1)/2;

disp('Starting MCMC.... ');
disp(' ' );
start_time = clock;
for loop = 1:nloop    
        %% sample beta
    invS = sparse(1:Tq,1:Tq,[1./Vbeta; repmat(invOmega22,T-1,1)]'); 
    K = H'*invS*H;   
    GinvOmega11 = bigG'*kron(speye(T),invOmega11);
    GinvOmega11G = GinvOmega11*bigG;
    invP = K + GinvOmega11G;
    C = chol(invP);                         % so that C'*C = invP
    betahat = C\(C'\(GinvOmega11*Y));
    beta = betahat + C\randn(Tq,1);         % C is upper triangular
        %% sample Omega11
    e1 = reshape(Y-bigG*beta,n,T);
    newS1 = S01 + e1*e1';
    invOmega11 = wishrnd(newS1\speye(n),newnu1);
    Omega11 = invOmega11\speye(n);
        %% sample Omega22
    e2 = reshape(H*beta,q,T);
    newS2 = S02 + sum(e2(:,2:end).^2,2)/2;
    invOmega22 = gamrnd(newnu2,1./newS2);
    Omega22 = 1./invOmega22;  
    if loop>burnin
        i=loop-burnin;
        store_beta(i,:) = beta';
        store_Omega11(i,:,:) = Omega11;
        store_Omega22(i,:) = Omega22';
    end    
    if ( mod( loop, 2000 ) ==0 )
        disp(  [ num2str( loop ) ' loops... ' ] )
    end    
end
disp( ['MCMC takes '  num2str( etime( clock, start_time) ) ' seconds' ] );
disp(' ' );

    %% compute posterior means
betahat = mean(store_beta)';
Omega11hat = squeeze(mean(store_Omega11)); 
Omega22hat = mean(store_Omega22); 
tempbetahat = reshape(betahat,nnp1,T)';
gamhat1 = tempbetahat(:,1:n+1);
gamhat2 = tempbetahat(:,(n+1)+1:2*(n+1));

The code gives no error when one lag is used, but when more than 1 lag is used, an error occurs reffering to the matrix "X" in line 28. The matrix X is basically constructed to include all the lags, in case of 6 lags the X matrix should have 12 columns with non-zero values, but when I use 6 lags, the Matrix "X" is having columns with zero values which generates the following error

Index in position 1 is invalid. Array indices must be positive integers or logical values.

The extended X matrix is genertated with the lines below, thus something has to be done here to avoid the error

for i=1:p
    X(:,(i-1)*n+1:i*n) = [Y0(3-i+1:end,:); shortY(1:T-i,:)];
end

But I don't know what should be done here, could anyone please help me with this issue?

[dpb]

Have you reached out to the authors of the code?

W/ no doc and no comments, makes it pretty hard-sledding for somebody to just look at the code and do anything in any reasonable amount of time...

[Ameer Fahmi]

I sent a message to the authors since few days and I haven't got any response, I couldn't wait more so I thought about sharing the problem here.

[Walter Roberson]

Y0 has 3 rows. p is maximum 6. for i=1:p means that i can be up to 6. Y0(3-i+1:end,:) when i reaches 4 gives 3-4+1 = 0 as a subscript.

Your p cannot exceed the number of rows in Y0. If you increase those rows then recode 3-i+1 as end-i+1

[Ameer Fahmi]

Is number of rows in Y0 is arbitrarily chosen? it can be any number porvided that the outcome of X matrix will not be zero, the authors in their original paper didn't collaborate on this point.

[Walter Roberson]

I do not know anything about the model.

Answer 1

[Geoff Hayes]

Ameer - if the lag p is 6, then in this code

for i=1:p
    X(:,(i-1)*n+1:i*n) = [Y0(3-i+1:end,:); shortY(1:T-i,:)];
end

when i is 4 (since iterating from 1 to 6), then

3 - i + 1 = 3 - 4 + 1 = 0

and so you are using an invalid index into Y0 and the error message makes sense. You need to adjust your code to handle this situation. Does Y0 need to be defined differently?

Y0 = USdata1(1:6,:);

And then adjust the 3 at

[Y0(3-i+1:end,:); shortY(1:T-i,:)];

?

[Ameer Fahmi]

Do you mean that the matrix Y0 is a matrix of the lags to be used? if we have 6 lags than the first 6 entries are taken in this matrix, I don't think this is the case, because the author of the code originally used one lag, and if you are right the matrix Y0 should be a vector of the first row, but the author used the three first entries in the Y0 matrix despite using one lag.

I got another error when I defined the Y0 matrix as you suggested.

[Geoff Hayes]

I see that there are only four columns in the USdata.csv spreadsheet, so what I suggested is not going to work. Looking at one of their papers, they are assuming models of 1 or 2 lags only...so having six is the problem because of

[Y0(3-i+1:end,:); shortY(1:T-i,:)];

Do you know why they are doing the above?

[Ameer Fahmi]

The X matrix is initially contsrtucted with 0 values for all the columns corresponding to the number of lags, the formula of X matrix with only zeros is referred to in the previous line

X = zeros(T,n*p);

Six lags and 2 variables will produce a matrix with 12 columns (6*2)

The formula you referred to is basically used to extend the matrix "ShortY" to include the number of lags and then pass it to the X matrix, "ShortY" is a matrix of the endogenous variables. However in the case of using 6 lags, the matrix "ShortY" which contains only the endogeous variables is not extended to encompass all the lags, which leaves the last columns with zero values when the matrix "ShortY" is passed to the matrix X, in my case the last six columns are still having zero values.

[dpb]

It appears you'll have to decipher the implementation of the algorithm and figure out what the original author should have done to handle more lags than his code does.

[Ameer Fahmi]

I'm trying to implemente the code line by line to figure out how it should be adjusted.

[dpb]

Good luck... :)