%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program estimates the generalized spatial balanced panel model put 
% forward by Baltagi, Egger and Pfaffermayr (2007) and Pfaffermayr (2009). 
% This is a test file and works with generated data. The matlab code is 
% based on Pfaffermayr(2009), Maximum Likelihood Estimation of a General 
% Unbalanced Spatial Random Effects Model - A Monte Carlo Study, Working 
% paper. This m-file refers to the unbalanced spatial panel.
%
% The follwing vectors and matrices have to be provided by the user.
% y...n x 1, dependent variable wit missings
% X...nx K, explanatory variabels including the constant
% So...n x NT, selector matrix for observed values
% Sm...m x NT, selector matrix for missings , m+n=N
% WN...NxN Spatial weighting matrix 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Functions in m-Files called in this procedure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPRE_MLE_unbal
% gendata for testing

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options and Organisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
warning off all;
minit=1;                 % Number valid ML-estimations
maxit=minit+2;           % Max number of ML-estimations
tol=1e-5;                % Tolerance 2 for two step likelihood maximization
options=optimset('Display','off',  'Diagnostics','off', 'MaxFunEvals',10000, ...
            'LargeScale','off', 'MaxIter',200,'OutputFcn', @outfun,...
            'TolFun',1e-5, 'TolX', 1e-5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters - to be skipped with true data  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta=[0.5;5];           % Slope and Constant
noind=size(beta,1);     % Number of slope parameters including the constant
vc=10;                  % Overall variance
theta=0.4;              % Share of cross sectional variance; s2m/s2n
s2n=vc*(1-theta);       % Remainder error
s2m=vc*theta;           % unit specific error
phi=s2m/s2n;            % Ratio of errorsd
rho1=0.8;               % Spatial correlation 1
rho2=0.8;               % spatial correlation 2
ss=0.3;                 % Percent missing
type=1;                 % 1...rook, 2... distance decay
cut=3;                  % Cut for distance based spatial weights
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Insert data here: y, X, So, Sm  WN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate data and weighing matrix for testing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=10^2;                  %Cross-section units
T=5;                    %Time dimension

% Generate data, the selection matrices and the spatial weighing matrix
[yt,Xt,WN,So,Sm] = gendata(N,T,beta,s2m,s2n,rho1,rho2,type,cut,ss);

% Variable names
vars=char('x1', 'cons');
K=size(vars,1);

% True model with missings
y=So*yt;
X=So*Xt;
clear  Xt yt 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Number of non-missings and missings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(So,1);
m=N*T-n;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ML-maximization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z0=[phi,rho1,rho2];   %starting values
tic
[estp,LogL,timer,iterv, V1md, V2md] = GSPRE_MLE_unbal(y,X,So,Sm,WN,N,T,z0,...
                                      minit,maxit,tol);
totaltime=toc;
fprintf('Total elapsed time %6.2f\n', totaltime)
fprintf('\n')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Standard errors of estimated paramters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sp1=diag(V1md).^0.5;
sp2=diag(V2md).^0.5;
sp=[sp1;sp2];
ttest=(estp(K+4,1)./sp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display ouput
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iter=mean(iterv);
ti=mean(timer);
time=ti/60;

% Variable names
par1=vars;
par2=char('s2m','s2n', 'rho1','rho2');
par3=char('Logl','n','T');

% Display computing time
disp('Iterations      min       sec')
disp([iter, time, ti])
fprintf('***************************************************************\n');

% Display estimation results
fprintf('\n');
fprintf('Unbalanced Panel\n');
fprintf('***************************************************************\n');
disp('         est. par.  standard error ' )
fprintf('***************************************************************\n');
for i=1:K
fprintf('%s %12.6f %12.6f\n',par1(i,:), estp(i,1),sp(i,1))  
end
fprintf('\n')
for i=1:4
fprintf('%s %12.6f %12.6f\n',par2(i,:), estp(K+i,1),sp(K+i,1))  
end
fprintf('\n')
for i=1:3
fprintf('%s %12.6f\n',par3(i,:), estp(K+4+i,1))  
end
fprintf('***************************************************************\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

