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

1 次查看（过去 30 天）
Ameer Fahmi
Ameer Fahmi 2019-7-18
评论： Walter Roberson 2019-7-19
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?
  5 个评论
显示 3更早的评论
Ameer Fahmi
Ameer Fahmi 2019-7-19
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Geoff Hayes
Geoff Hayes 2019-7-18
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,:)];
?
  6 个评论
显示 4更早的评论

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by