function [P,VC_U,VC_pv1,pvec,std_pvec,t_pvec,WaldStat,FStat]=sure(Y,X0,const,sigma)
% This function estimates the parameters of the SUR system
% as well as test statistics for the significance of the intercepts

% Inputs: Y                - matrix of the dependent variables, 
%            X0               - matrix of explanatory variables 
%            const={0,1}  - if chosen 1 includes the intercept, if 0 -no intercept
%            sigma          - true variance covariance matrix 
% Outputs: P              -  estimated coefficients, 
%               VC_U       - variance-covariance matrix of the residuals, 
%               VC_pv      - variane-covariance matrix for the estimated coefficients, 
%               pvec         - vectorized estimates of the coefficients, 
%               std_pvec   - vectorized standard deviations of the estimated coefficients, 
%               t_pvec       - vectorized values of the t-statistic for the significance of the coefficients,
%               WaldStat   - value of the test statistic for a Wald-Test of
%                                 joint hypothesis of significance of the
%                                 intercept, using the true
%                                 variance-covariance matrix sigma
%               FStat         - value of the test statistic for the F-Test of
%                                 joint hypothesis of significance of the
%                                 intercept, using the estimated VC matrix
%                                
% The notation follows Simmet and Pohlmeier (2018)
% This version 02/2019

[t k]=size(X0);      % returns the dimensions of X0
[t n]=size(Y);       % returns the dimensions of Y
if const==1
  X=[ones(t,1) X0];   % adds the column of ones if intercept is required
  k=k+1;
else 
 X=X0; 
 k=k;
end
  P=inv(X'*X)*X'*Y;       % estimated coefficients
  Y_hat = X*P;           % fitted values for Y
  u = Y - Y_hat;          % residuals
  VC_U = u'*u/(t-k); 
  VC_pv1 =kron(VC_U,inv(X'*X));
     
 % t-statistics for the test of the significance of the coefficients    
  pvec = reshape(P,[],1);
  ttest = reshape(P-[zeros(1,n); ones(1,n)],[],1);
  std_pvec = sqrt(diag(VC_pv1));
  t_pvec = (ttest)./ std_pvec;   
  % F- test statistics for joint significance of the intercepts
   R_alpha=zeros(n,2*n);
for i=1:n 
   R_alpha(i,2*i-1)=1; 
end;
   alpha=R_alpha*pvec;
   VC_alpha1 = R_alpha * VC_pv1 * R_alpha';
   FStat = ((t-n-1)/n/t)*alpha'*inv(VC_alpha1) * alpha;

  % Wald-test statistic 
   VC_pv2 =kron(sigma,inv(X'*X));
   VC_alpha2 = R_alpha * VC_pv2 * R_alpha';
   WaldStat = alpha'*inv(VC_alpha2) * alpha;
end
