MATRIX ARRANGEMENT PROBLEM MATLAB

Question

Saptorshee Chakraborty 2022-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1740245-matrix-arrangement-problem-matlab

编辑： Cris LaPierre 2022-6-16

在 MATLAB Online 中打开

Hello,

I am working with panel structural breaks, where , and X1, X2, X3 (d=3).

Now I have three breaks, B=3, while the author uses the code with B=1 and GB, GA (for me it will be GI, GII, GIII, GIV)

I am stuck at T_frac, Y_frac, X_frac somehow the index becomes more than 3.

I am permission to use the code of the author, any help will be appreciated.

Xvec                    = [x1, x2, x3];
Yvec                    = y;
N=114;
TT=48;
xlag                    = x(:,1:TT-1,:);
ylag                    = y(:,2:TT);
T                       = TT-1;
Xmat                    = xlag;
ymat                    = ylag;
p       = size(Xmat,3);
S       = 30;                                                               
SSR     = zeros(1,(T-2));
SSR(1)  = inf;
ymatfd  = diff(ymat,1,2);
Xmatfd  = diff(Xmat,1,2);
for j = 1:B
    for bk = 2:(T-2)
        regime = [0,bk,T-1];                                                % Get the bounds of potential subperiods      
        for aa = 1:(j+1)                             
            T_frac          = length((regime(aa)+1):(regime(aa+1)));            
            y_frac          = reshape(ymatfd(:,(regime(aa)+1):(regime(aa+1)))',N*T_frac,1);            
            for l = 1:p            
                X_frac(:,l) = reshape(Xmatfd(:,(regime(aa)+1):(regime(aa+1)),l)',N*T_frac,1);                
            end            
            time_frac       = repmat(1:T_frac,1,N)';            
            code_frac       = kron((1:N)',ones(T_frac,1));                  
            idxdata_frac    = dataset(code_frac, time_frac, y_frac, X_frac);                        
            idxdata_frac.Properties.VarNames = {'N'  'T'  'y'  'X'};
            if aa == 1
                [~,~,ssr,]  = est_group(idxdata_frac,GB);
            elseif aa == 2
                [~,~,ssr,]  = est_group(idxdata_frac,GA);
            end                
            SSR(bk)         = SSR(bk) + ssr;        
            clear X_frac y_frac idxdata_frac time_frac code_frac
        end
    end
    [resQ,k]       = min(SSR);                                              % Compare SSR to get a breakdate
end

2 个评论
显示无隐藏无

Jan 2022-6-14

"Now I have three breaks, B=3, while the author uses the code with B=1 and GB, GA (for me it will be GI, GII, GIII, GIV)" - remember, that the readers do not know, what you are talking off. Who is the autheor? What us GB and GA, GI, GII, GIII and GIV?

What does this mean: "I am stuck at T_frac, Y_frac, X_frac somehow the index becomes more than 3."?

Saptorshee Chakraborty 2022-6-14

When I compute using B=3.

请先登录，再进行评论。

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

Answer 1

Cris LaPierre 2022-6-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1740245-matrix-arrangement-problem-matlab#answer_985835

在 MATLAB Online 中打开

Not that your loop counter is j, but your indexing uses j+1, and then again use loop counter aa, but index using aa+1. Therefore, when B=3, j+1=4 and aa+1=5.

for j = 1:B
    ...
        for aa = 1:(j+1)
            T_frac          = length((regime(aa)+1):(regime(aa+1)));
            %                                               ^^^^
        end
end

12 个评论
显示 10更早的评论隐藏 10更早的评论

Saptorshee Chakraborty 2022-6-15

在 MATLAB Online 中打开

The Author's code is of three pages, I am pasting it on a single page. If I use B=1 with my data, the CODE works fine, but if I use with B=3 and change GB and GA to GI=3, GII=3, GIII=2, GIV=3 and changes in page 1,2 and 3 the code fails with the erroIndex exceeds the number of array elements. Index must not exceed 3.

Error in est_gpbk_fd_tvG (line 39)

T_frac = length((regime(aa)+1):(regime(aa+1)));

Error in application (line 43)

[k_lsbg,gpI_lsbg,gpII_lsbg,gpIII_lsbg,gpIV_lsbg,bI_lsbg,bII_lsbg,bIII_lsbg,bIV_lsbg,seI_lsbg,seII_lsbg,seIII_lsbg,seIV_lsbg,~] = est_gpbk_fd_tvG(ymat,Xmat,B,GI,GII,GIII,GIV);

Any help is appreciated

%Page 1
load('data.mat');
Xvec                    = [X1 X2 X3];                   
Yvec                    = y;
TT                      = 48;
N                       = length(Yvec)/TT;
K                       = size(Xvec,2);
for k = 1:K
    x(:,:,k)            = reshape(Xvec(:,k), TT, N)';
end
y                       = reshape(Yvec,TT,N)';
xlag                    = x(:,1:TT-1,:);
ylag                    = y(:,2:TT);
T                       = TT-1;
Xmat                    = xlag;
ymat                    = ylag;
GB                      = 5;                                                
GA                      = 7;                                                
B                       = 1; 
[k_lsbg,gpB_lsbg,gpA_lsbg,bB_lsbg,bA_lsbg,seB_lsbg,seA_lsbg,~]   = est_gpbk_fd_tvG(ymat,Xmat,B,GB,GA);
%PAGE- 2
function[k,gpmemB,gpmemA,coefB,coefA,seB,seA,resQ] = est_gpbk_fd_tvG(ymat,Xmat,B,GB,GA)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This method simultaneously estimate the group structure, break point, and
% slope coefficients for panel models with fixed effects. 
% It allows that the number of groups varies over time.
%
% Input: 
% ymat: dependent variable in matrix form
% Xmat: explanatory variables in matrix form
% B: number of breaking points
%
% Output:
% k: estimated break point
% gpmemB: group membership estimate before the break
% gpmemA: group membership estimate after the break
% coefB: coefficient estimates before the break
% coefA: coefficient estimates after the break
% seB: standard error before the break
% seA: standard error after the break
% resQ: sum of squared residuals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T       = size(ymat,2);
N       = size(ymat,1);
p       = size(Xmat,3);
S       = 30;                                                              
SSR     = zeros(1,(T-2));
SSR(1)  = inf;
ymatfd  = diff(ymat,1,2);
Xmatfd  = diff(Xmat,1,2);
%% Break point estimation
for j = 1:B
    for bk = 2:(T-2)
        regime = [0,bk,T-1];                                                % Get the bounds of potential subperiods      
        for aa = 1:(d+1)   ;                             
            T_frac          = length((regime(aa)+1):(regime(aa+1)));            
            y_frac          = reshape(ymatfd(:,(regime(aa)+1):(regime(aa+1)))',N*T_frac,1);            
            for l = 1:p            
                X_frac(:,l) = reshape(Xmatfd(:,(regime(aa)+1):(regime(aa+1)),l)',N*T_frac,1);                
            end            
            time_frac       = repmat(1:T_frac,1,N)';            
            code_frac       = kron((1:N)',ones(T_frac,1));                  
            idxdata_frac    = dataset(code_frac, time_frac, y_frac, X_frac);                        
            idxdata_frac.Properties.VarNames = {'N'  'T'  'y'  'X'};
            if aa == 1
                [~,~,ssr,]  = est_group(idxdata_frac,GB);
            elseif aa == 2
                [~,~,ssr,]  = est_group(idxdata_frac,GA);
            end                
            SSR(bk)         = SSR(bk) + ssr;        
            clear X_frac y_frac idxdata_frac time_frac code_frac
        end
    end
    [resQ,k]       = min(SSR);                                              % Compare SSR to get a breakdate
end
regime          = [0,k,T-1];
coefB           = zeros(GB,p);
coefA           = zeros(GA,p);
seB             = zeros(GB,p);
seA             = zeros(GA,p);
gpmem           = [];
%% Post break estimation: Classification for each subperiod
for j = 1:(B+1)            
    for aa = 1:j                                 
        T_frac          = length((regime(aa)+1):(regime(aa+1)));                    
        y_frac          = reshape(ymatfd(:,(regime(aa)+1):(regime(aa+1)))',N*T_frac,1);                   
        for l = 1:p                    
            X_frac(:,l) = reshape(Xmatfd(:,(regime(aa)+1):(regime(aa+1)),l)',N*T_frac,1);                            
        end        
        time_frac       = repmat(1:T_frac,1,N)';        
        code_frac       = kron((1:N)',ones(T_frac,1));                          
        idxdata_frac    = dataset(code_frac, time_frac, y_frac, X_frac);                                
        idxdata_frac.Properties.VarNames = {'N'  'T'  'y'  'X'};
        if aa == 1
            [coefB,gpmemB,~,seB]     = est_group(idxdata_frac,GB);
        elseif aa == 2
            [coefA,gpmemA,~,seA]     = est_group(idxdata_frac,GA);
        end
        clear X_frac y_frac idxdata_frac time_frac code_frac
    end
end
%Page 3
function [est_beta, est_group, resQ_best, est_se] = est_group(idxdata,G)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function estimates the group memberships and slope coefficients
% for panel models without structural breaks
%
% Input:
% idxdata: data set
% G: number of groups
%
% Output:
% est_beta: estimated slope coefficients
% est_group: estimated group memberships
% resQ_best: sum of squared residuals
% est_se: estimated standard error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N                               = idxdata.N(end); 
T                               = size(idxdata,1)/N;
K                               = size(idxdata.X,2);
Y                               = idxdata.y;
X                               = idxdata.X;
%% initialization
gi_auxaux                       = zeros(N,G);
if G == 1
    Nsim                        = 1;                                        
else 
    Nsim                        = 1000;                                     % number of initial values (specified by the authors)
end
resQ_best                       = realmax;                                  % sum of squared residuals for the best initial values
max_iteration                   = 20;
%% estimation
for jsim = 1:Nsim                                                           % use different initial values 
    gi_init                     = zeros(N,G);                                                   
    for i = 1:N                                                             % initialize a random group pattern
        perm                    = randperm(G);    
        for g = 1:G    
            gi_init(i,g)        = (perm(1)==g);        
        end    
    end
    gi                          = gi_init;    
    par_init                    = zeros(N*T,K);
        
    deltapar                    = 1;                                        % iteration indicator   
    num_iter                    = 0;                                        % number of iterations
    while deltapar > 0 && num_iter < max_iteration                          % iterated optimization
        if ~all(sum(gi))
            for i = 1:N                                                     % if any group has no element, start a new random initialization
            perm                = randperm(G);    
                for g = 1:G    
                    gi_init(i,g) = (perm(1)==g);        
                end    
            end
            gi                  = gi_init;           
        end
        %% step 1: update          
        for g = 1:G 
            NN                  = 1:N;
            gi                  = logical(gi);
            this_group          = gi(:,g);    
            g_index             = NN(this_group);    
            g_data              = idxdata(ismember(idxdata.N,g_index),:);   % group-specific data	
            Ng                  = sum(gi(:,g));                             % number of individuals within this group
            gX                  = g_data.X;
            gy                  = g_data.y;
            alpha               = (gX'*gX)\(gX'*gy);                        % regression for each group
            beta(g,:)           = alpha';
            res                 = gy-gX*alpha;
            xi                  = repmat(res,1,K).*gX;   
            Omega               = (xi'*xi)/(Ng*T);
            Phi                 = (gX'*gX)/(Ng*T)\eye(K);
            S                   = Phi*Omega*Phi'/(T*Ng);
            se(g,:)             = reshape(sqrt(diag(S)),1,K);                        
        end
        
        %% step 2: group assignment
        for g = 1:G    
            beta_group_temp     = beta(g,:);        
            beta_group          = kron(ones(N*T,1),beta_group_temp);        % beta of the g-th group     
            U                   = (Y-sum(X.*beta_group,2)).^2;              % squared residual using beta of the g-th group  
            RU                  = reshape(U,T,N)';                   
            gi_auxaux(:,g)      = sum(RU')';                                % time summation of squared residual for the g-th group
        end    
        [gi_class, gindt]       = min(gi_auxaux,[],2);                      % for each individual, find the group with least sum of squared residual
        for g = 1:G
            gi(:,g)             = (gi_auxaux(:,g) == gi_class);             % assign each individual to the optimal group       
        end
        %% check convergence
        par_new                 = zeros(N*T,K);                       
        for i = 1:N
            par_new_indiv       = beta(gindt(i),:);                         % updated parameter for each individual 
            par_new((i-1)*T+1:i*T,:) = repmat(par_new_indiv,T,1);           % updated parameter
        end
        deltapar                = norm(par_new-par_init);                   % update iteration indicator: check norm between updated and initial parameter
        par_init                = par_new;    
        num_iter                = num_iter + 1;
    end
    resQ                        = (Y-sum(X.*par_new,2))'*(Y-sum(X.*par_new,2));                  
    if resQ < resQ_best
        resQ_best               = resQ;
        gi_best                 = gi;
        beta_best               = beta;
        se_best                 = se;
    end        
    jsim
end
%% final estimate
est_group                       = gi_best;
est_beta                        = beta_best;
est_se                          = se_best;
para                            = par_new;
end

Cris LaPierre 2022-6-16

编辑：Cris LaPierre 2022-6-16

I didn't mention a fix. I just pointed out the code that was causing the problem.

The fix is to make regime longer. It needs to be a 1x(B+2) vector to not error. Again, it is currently 1x3, so when B=1, the code works. When B=3, it needs to be 1x5.

I don't know what values to add, though. That is something where you would need to apply your knowledge in this space to corectly adjust the code.

Saptorshee Chakraborty 2022-6-16