Asked by Ameer Fahmi
on 18 Jul 2019

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?

Answer by Geoff Hayes
on 18 Jul 2019

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,:)];

?

dpb
on 18 Jul 2019

Ameer Fahmi
on 19 Jul 2019

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

dpb
on 19 Jul 2019

Good luck... :)

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 5 Comments

## dpb (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/472343-error-index-in-position-1-is-invalid-array-indices-must-be-positive-integers-or-logical-values#comment_725748

## Ameer Fahmi (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/472343-error-index-in-position-1-is-invalid-array-indices-must-be-positive-integers-or-logical-values#comment_725809

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/472343-error-index-in-position-1-is-invalid-array-indices-must-be-positive-integers-or-logical-values#comment_725940

## Ameer Fahmi (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/472343-error-index-in-position-1-is-invalid-array-indices-must-be-positive-integers-or-logical-values#comment_726072

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/472343-error-index-in-position-1-is-invalid-array-indices-must-be-positive-integers-or-logical-values#comment_726335

Sign in to comment.