% Function to compute point-optimal statistic from Moon, Perron, Phillips (2007) with fixed initial condition

% Usage: MPP = mpp(z,p,cset)

% where    z: matrix of panel data of dimension T x N
%          p: deterministic component (p=-1 for no deterministic component, p=0 for constant and p = 1 for constant and linear trend)
%          cset: vector of length N of local to unity parameters for computing statistic

% Date: April 1st, 2005

function mppstat=mpp(z,p,cset);

[t,n] = size(z);

                        
                        %compute likelihood under null and alternative

                        lnc = zeros(2,1);
                        
                        for cindex = 1:2

                            if cindex == 1
                                ci = zeros(n,1);
                            else
                                ci = cset;   
                            end
                            

                           % quasi-differences

                                if p == 1;
                                    rhoc = ones(n,1)-ci/(n^0.25*t);
                                    g = [ones(t,1),[1:t]'];
                                    btilde = zeros(n,2);                        
                               
                                else
                                    rhoc = ones(n,1)-ci/(n^0.5*t);
                                    g = ones(t,1);
                                    btilde = zeros(n,1);                           
                                end

                                for j = 1:n
                                    deltaz(:,j) = [z(1,j);z(2:t,j)-z(1:t-1,j)*rhoc(j)];                                
                                end
                                ztilde=zeros(t,n);
                           
                           if p >=0;
                                %remove deterministic components by GLS 
                            
                                for j=1:n
                                
                                    deltag = [g(1,:);g(2:t,:)-g(1:t-1,:)*rhoc(j)];
                                    btilde(j,:) = (inv(deltag'*deltag)*deltag'*deltaz(:,j))';
                                    ztilde(:,j) = deltaz(:,j) - deltag*btilde(j,:)';
                                    
                                end
                                
                            else
                               ztilde = deltaz;
                            end
                            
                        %estimate sigma
                        
                        if cindex == 1;
                            sigmahat =  sum(ztilde.^2)/t;
                        end
                        
                         % compute sum of squared residuals 
                                                    
                            Ztilde = reshape(ztilde./(ones(t,1)*sqrt(sigmahat)),t*n,1);

                            lnc(cindex) = Ztilde'*Ztilde;
                            
                        end
                     
                        % compute centering (omega) terms if trend
                        
                        if p ==1;
                            trend2 = [0:1/t:1-1/t]';
                            omega2 = -sum(trend2)/t + 2/t*trend2'*trend2 - 1/3;
                        
                            sum1 = 0;                       
                            for t1=1:t
                                terms = (trend2(t1).*trend2).*min(trend2(t1),trend2);
                                sum1= sum1 + sum(terms);
                            end
                        
                            omega4 = sum1/t^2 - 2/(3*t)*trend2'*trend2 + 1/9;
                        end
    
                        % compute statistic
                        
                        if p == 1
                            mppstat = sqrt(45/mean(cset.^4))*(lnc(2)-lnc(1) +sum(cset)/n^0.25+sum(cset.^2)/sqrt(n)*omega2+sum(cset.^4)/n*omega4);
                        else
                            mppstat = sqrt(0.5/mean(cset.^2))*(lnc(2)-lnc(1) -0.5*mean(cset.^2));
                        end