% Code to estimate Linear Expenditure System with endogenous interaction
% effects and cost differences among jurisdictions
%
% written by: J.Paul Elhorst Winter 2010/2011
% University of Groningen
% Department of Economics
% 9700AV Groningen
% the Netherlands
% j.p.elhorst@rug.nl
%
% REFERENCE:
% Allers M.A., Elhorst J.P. (2011) A simultaneous equations model of fiscal 
% policy interactions. Journal of Regional Science. Forthcoming.
%
% Read data and spatial weights matrix
D=wk1read('y:\data\data.wk1',1,0);
W=wk1read('y:\data\BC.wk1');
N=496;
nobs=N;
W=normw(W); % row-normalize W, routine LeSage
%
% Specify the Linear Expenditure system. The specification below is used in
% the Journal of Regional Science paper. Since your model will be
% different, the specification below needs to be adjusted. It is
% recommended to start with a simple model first, e.g., 3 equations, and N
% is small, e.g., N<100. If the program works and is able to find a
% maximum, then the model can be extended
% _________________________________________________________________________
neq=7;                      % number of equations in the system
% NOTE: Many lines below NEED TO BE ADJUSTED IF neq IS NOT EQUAL TO 7;
% these lines are indicated
m=neq-1;                    % Note one equation is eliminated because of 
% the adding-up restrictions (see explanation below eq. 22). This program
% eliminates the last equation.

y=zeros(nobs,neq);

y(:,1)=D(:,76); % tax revenue OZB
y(:,2)=D(:,64); % welfare BIJSTAND
y(:,3)=D(:,65); % education, health & social services ZORG
y(:,4)=D(:,66)+D(:,67); %Arts, leisure, green space & open air recreation ALGO
y(:,5)=D(:,70)+D(:,71); %Fire brigade, civil protection & environment FCE
y(:,6)=D(:,75)+D(:,77); %Roads, waterways, monuments & historical objects RWMH
y(:,7)=D(:,72); %Administration & register office REST

y(:,1)=-y(:,1);             % negative in case of revenues instead of expenditures
xi=sum(y,2);                % sum of all expenditures minus taxes

% After the %-sign, you can find the (Dutch) names of the explanatory
% variables we used in this study

xbeta=D(:,[79,22,24,40]);   % constant, rechts, inkhh en taxpr
x1= D(:,[79,32,27]);        % OZB: constant, grondslag, lihh
x2= D(:,[79,28,38]);        % BIJSTAND: constant, bijstontv, uitkossd
x3= D(:,[79,35,30,39,34]);  % ZORG: constant, leerl, kpreg, uitkontv, minderh
x4= D(:,[79,33,29,30]);     % ALGO: constant, woonruim, kplok, kpreg
x5= D(:,[79,31,47,36]);     % FCE:  constant, oad, bedrvest, oppbeb
x6= D(:,[79,54,36,45]);     % RWMH: constant, oeverlbodem, oppbeb, opphistk 
x7=D(:,[79,30,33,34]);      % REST: constant, kpreg, woonruim, minderh
x=[x1 x2 x3 x4 x5 x6 x7];

wy=W*y;
wx=W*x;

npar= [cols(x1);cols(x2);cols(x3);cols(x4);cols(x5);cols(x6);cols(x7)]; % number of gam parameters in each equation
cumpar=[0;npar(1);npar(1)+npar(2);npar(1)+npar(2)+npar(3);npar(1)+npar(2)+npar(3)+npar(4);npar(1)+npar(2)+npar(3)+npar(4)+npar(5);npar(1)+npar(2)+npar(3)+npar(4)+npar(5)+npar(6);npar(1)+npar(2)+npar(3)+npar(4)+npar(5)+npar(6)+npar(7)]; % number of gam parameters in npar cumulated

% Names of the variables used when printing the results

vLES=strvcat('LES');
v1beta=strvcat('OZB-beta',     'OZB-rechts',     'OZB-inkhh',     'OZB-taxpr');
v2beta=strvcat('BIJSTAND-beta','BIJSTAND-rechts','BIJSTAND-inkhh','BIJSTAND-taxpr');
v3beta=strvcat('ZORG-beta',    'ZORG-rechts',    'ZORG-inkhh',    'ZORG-taxpr');
v4beta=strvcat('ALGO-beta',   'ALGO-rechts',   'ALGO-inkhh',   'ALGO-taxpr');
v5beta=strvcat('FCE-beta',     'FCE-rechts',     'FCE-inkhh',     'FCE-taxpr');
v6beta=strvcat('RWMH-beta',    'RWMH-rechts',    'RWMH-inkhh',    'RWMH-taxpr');
v7beta=strvcat('REST-const','REST-rechts', 'REST-inkhh', 'REST-taxpr');

v1names=strvcat('OZB-const','OZB-grondsl','OZB-lihh');
v2names=strvcat('BIJSTAND-const','BIJSTAND-bijstontv','BIJSTAND-uitkossd');
v3names=strvcat('ZORG-const','ZORG-leerl','ZORG-kpreg','ZORG-uitkontv','ZORG-minderh');
v4names=strvcat('ALGO-const','ALGO-woonruim','ALGO-kplok','ALGO-kpreg');
v5names=strvcat('FCE-const','FCE-oad','FCE-bedrvest','FCE-oppbeb');
v6names=strvcat('RWMH-const','RWMH-oeverlbodem','RWMH-oppbeb','RWMH-opphistk');
v7names=strvcat('REST-const','REST-kpreg','REST-woonruim','REST-minderh');

dnames=strvcat('delta1','delta2','delta3','delta4','delta5','delta6','delta7');

% end of specification of the linear expenditure system
% _________________________________________________________________________
info.lflag=0;
[rmin,rmax,convg,maxit,detval,ldetflag,eflag,order,miter,options] = sem_parse(info); % routine LeSage
[rmin,rmax,time2] = sem_eigs(eflag,W,rmin,rmax,N); % routine LeSage
[detval,time1] = sem_lndet(ldetflag,W,rmin,rmax,detval,order,miter); % routine LeSage
fid=1;
[s1 s2]=size(xbeta);
nobeq=s2;
nobeta=m*nobeq;
beta=zeros(nobeta,1);
for i=1:nobs
    xib(i,:)=xi(i)*xbeta(i,:); %interaction effects between variables X explaining discretionary income and variables S explaining spending needs
end
[s1 s2]=size(x);
nogam=s2;
gam=zeros(nogam,1);
par=[beta;gam];
totpar=length(par);

delta=zeros(neq,1);
% cross-product matrix of xbeta and x=[x1 x2 x3]
xz=zeros(nobs,nobeq*nogam);
wxz=zeros(nobs,nobeq*nogam);
for i=1:nobs
for k=1:nogam
    for j=1:nobeq
        xz(i,(k-1)*nobeq+j)=xbeta(i,j)*x(i,k);
        wxz(i,(k-1)*nobeq+j)=xbeta(i,j)*wx(i,k);
    end;
end;
end;
emat = zeros(nobs,neq);
xc=ones(N,1);
index=zeros(nogam,1);
for k=1:nogam
    for j=1:neq
        if (k>cumpar(j) & k < cumpar(j+1)+1 ) index(k)=j; end
    end
end
% ************************OLS of each single equation *********************
for i=1:neq
    if (i<neq) 
        results=ols(y(:,i),[xib x(:,cumpar(i)+1:cumpar(i+1))]);
        if (i==1) vnames=strvcat(vLES,v1beta,v1names);
        elseif (i==2) vnames=strvcat(vLES,v2beta,v2names);
        elseif (i==3) vnames=strvcat(vLES,v3beta,v3names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        elseif (i==4) vnames=strvcat(vLES,v4beta,v4names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        elseif (i==5) vnames=strvcat(vLES,v5beta,v5names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        else vnames=strvcat(vLES,v6beta,v6names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        end
        prt_reg(results,vnames);
        par((i-1)*nobeq+1:i*nobeq)=results.beta(1:nobeq);
        par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1))=results.beta(nobeq+1:nobeq+npar(i));
        emat(:,i)=results.resid;
        sigma(i,i)=results.si2;
    else 
        results=ols(y(:,i),[xib x(:,cumpar(i)+1:cumpar(i+1))]);
        vnames=strvcat(vLES,v7beta,v7names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        prt_reg(results,vnames);
        par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1))=results.beta(nobeq+1:nobeq+npar(i));
        emat(:,i)=results.resid;
        sigma(i,i)=results.si2;
    end;
end;
sigmavec=zeros(0.5*m*(m+1),1);
tel=0;
for i=1:m
    for j=i:m
        tel=tel+1;
        sigmavec(tel)=sigma(i,j);
    end
end
ML_ols=-f2_lessarquotient([par;zeros(neq,1);sigmavec],nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,W)
% ************************ML estimation of each equation including spatially lagged dependent variable *********************
%
% sardetval is rewritten version of LeSage's routine "sar" and is meant to
% save computation time and to work with one constant matrix "detval" of
% the Jacobian term in the log-likelihood function throughout the whole
% computer run
%
for i=1:neq
    if (i<neq) 
        results=sardetval(y(:,i),[xib x(:,cumpar(i)+1:cumpar(i+1))],W,info,detval,rmin,rmax,ldetflag,miter,order);
        if (i==1) vnames=strvcat(vLES,v1beta,v1names);
        elseif (i==2) vnames=strvcat(vLES,v2beta,v2names);
        elseif (i==3) vnames=strvcat(vLES,v3beta,v3names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        elseif (i==4) vnames=strvcat(vLES,v4beta,v4names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        elseif (i==5) vnames=strvcat(vLES,v5beta,v5names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        else vnames=strvcat(vLES,v6beta,v6names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        end
        prt_spat(results,vnames);
        deltasar(i,1)=results.rho;
        parsar((i-1)*nobeq+1:i*nobeq)=results.beta(1:nobeq);
        parsar(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1))=results.beta(nobeq+1:nobeq+npar(i));
        ematsar(:,i)=results.resid;
        sigmasar(i,i)=results.sige;
    else 
        results=sardetval(y(:,i),[xib x(:,cumpar(i)+1:cumpar(i+1))],W,info,detval,rmin,rmax,ldetflag,miter,order);
        vnames=strvcat(vLES,v7beta,v7names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
        prt_spat(results,vnames);
        deltasar(i,1)=results.rho;
        parsar((i-1)*nobeq+1:i*nobeq)=results.beta(1:nobeq);
        parsar(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1))=results.beta(nobeq+1:nobeq+npar(i));
        ematsar(:,i)=results.resid;
        sigmasar(i,i)=results.sige;
    end;
end;
sigmavec=zeros(0.5*m*(m+1),1);
tel=0;
for i=1:m
    for j=i:m
        tel=tel+1;
        sigmavec(tel)=sigmasar(i,j);
    end
end
ML_sar=-f2_lessarquotient([parsar';deltasar;sigmavec],nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,W)
sumpar=zeros(nobeq,1);
for j=1:m
    j1=(j-1)*nobeq+1;
    j2=j*nobeq;
    sumpar=sumpar+parsar(j1:j2)';
end
betafactor=zeros(nobs,neq);
for j=1:neq-1
    j1=(j-1)*nobeq+1;
    j2=j*nobeq;
    for i=1:nobs
        betafactor(i,j)=xbeta(i,:)*parsar((j-1)*nobeq+1:j*nobeq)';
    end
end
j=neq;
for i=1:nobs
    betafactor(i,j)=xbeta(i,:)*sumpar(1:nobeq);
end
% ***********************Estimation of LES without spatial effects ******
itermax=100;
vnames=strvcat(vLES,v1beta,v2beta,v3beta,v4beta,v5beta,v6beta,v1names,v2names,v3names,v4names,v5names,v6names,v7names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
% initial sigma matrix using ols
sigma = zeros(m,m);
for i=1:m
 for j=i:m
    sigma(i,j) = (emat(:,i)'*emat(:,j))/nobs;
    if j > i; sigma(j,i) = sigma(i,j); end
 end;
end;
sigmai=inv(sigma)
ypar=zeros(m,1,nobs); 
xpar=zeros(m,totpar,nobs); % m is number of equations in system (one eq. eliminated), totpar is number of parameters, nobs is number of observations
options.Display='off';
options.MaxFunEvals=1000;
options.MaxIter=1000;
options.TolX=0.001;
options.TolFun=0.001;
iter=0;
converge=1.0;
criteria=0.001;
parit=zeros(totpar,itermax);

for i=1:nobs % ypar en xpar with respect to beta parameters
    for j=1:m
        ypar(j,1,i)=y(i,j);
            for k=1:nobeq
                p=(j-1)*nobeq+k;
                xpar(j,p,i)=xib(i,k);
            end
    end;
end;

while ( converge>criteria & iter<itermax)
iter=iter+1;
parsigold=par;
%
% Alternately estimate gamma and beta parameters
%
iter1=0;    
converge1=1.0;
while (converge1>criteria & iter1<itermax)
iter1=iter1+1;
sumpar=zeros(nobeq,1);
for j=1:m
    j1=(j-1)*nobeq+1;
    j2=j*nobeq;
    sumpar=sumpar+par(j1:j2);
end
for i=1:nobs
    for j=1:m
        j1=(j-1)*nobeq+1;
        j2=j*nobeq;
        for k=1:nogam
            p=nobeta+k;
            k1=(k-1)*nobeq+1;
            k2=k*nobeq;
            if (k > cumpar(j) & k < cumpar(j+1)+1 ) 
                xpar(j,p,i)=x(i,k)-xz(i,k1:k2)*par(j1:j2);
            else xpar(j,p,i)=-xz(i,k1:k2)*par(j1:j2);
            end;
        end;
    end;
end;
xvx=zeros(totpar,totpar);
xvy=zeros(totpar,1);
for i=1:nobs
    yhulp(:,:)=ypar(:,:,i);
    xhulp(:,:)=xpar(:,:,i);
    xvx=xvx+xhulp'*sigmai*xhulp;
    xvy=xvy+xhulp'*sigmai*yhulp;
end;
parold=par;
par=inv(xvx)*xvy;
converge1=sum(abs(par-parold));
end; % end loop alternately estimate gamma and beta parameters

itertot(iter)=iter1;
parit(:,iter)=par;
converge=sum(abs(par-parsigold));
for i=1:nobs
    eres=ypar(:,:,i)-xpar(:,:,i)*par;
    emat(i,1:m)=eres';
end
for i=1:m
 for j=i:m
    sigma(i,j) = (emat(:,i)'*emat(:,j))/nobs;
    if j > i;sigma(j,i) = sigma(i,j);end
 end;
end;
sigmai=inv(sigma);
end;
iter
itertot
%
% Measures of R-squared
%
for i=1:nobs
    ymat(i,:)=y(i,:)-mean(y);
end
xt=emat'*emat;
xn=ymat(:,1:m)'*ymat(:,1:m);
xtel=0;
xnoem=0;
for i=1:nobs
    xtel=xtel+emat(i,1:m)*inv(sigma(1:m,1:m))*emat(i,1:m)';
    xnoem=xnoem+ymat(i,1:m)*inv(sigma(1:m,1:m))*ymat(i,1:m)';
end
for i=1:m
    R2(i)=1-xt(i,i)/xn(i,i);
end
R2sys1=1-xtel/xnoem;
xtr=diag(xt);
R2sys2=R2*xtr(1:m)/sum(xtr(1:m));
%
% Calculation var-cov matrix for T-values
%
pp=0.5*m*(m+1);
xpx=zeros(totpar,totpar);
sigmami=inv(sigma(1:m,1:m));
% par, par (beta and gamma)
tpar=zeros(m,totpar,nobs);
for i=1:nobs
    for j=1:m
        for k=1:nobeq
            p=0;
            for t1=1:nogam
                t2=k+(t1-1)*nobeq;
                p=p+xz(i,t2)*par(nobeta+t1);
            end
            pk=(j-1)*nobeq+k;
            tpar(j,pk,i)=xib(i,k)-p;
        end
        for k=nobeta+1:totpar
            tpar(j,k,i)=xpar(j,k,i);
        end
    end
end
xtx=zeros(totpar,totpar);
for i=1:nobs
    xhelp(:,:)=tpar(:,:,i);
    xtx=xtx+xhelp'*sigmami*xhelp;
end
xpx(1:totpar,1:totpar)=xtx;
xpxi=inv(xpx);
xdiag=diag(xpxi);
bout=par;
tstat=bout./sqrt(xdiag(1:totpar));
%print results
fid=1;
fprintf(fid,'\n');
for i=1:m
    fprintf(fid,'R-squared eq.   = %2d,%9.4f   \n',i,R2(i));
end
fprintf(fid,'R-squared system McElroy en Dhrymes  = %11.4f,%11.4f  \n',R2sys1,R2sys2);
% now print coefficient estimates, t-statistics and probabilities
tout = norm_prb(tstat); % find asymptotic z (normal) probabilities
tmp = [bout tstat tout];  % matrix to be printed
% column labels for printing results
bstring = 'Coefficient'; tstring = 'Asymptot t-stat'; pstring = 'z-probability';
cnames = strvcat(bstring,tstring,pstring);
in.cnames = cnames;
in.rnames = vnames;
in.fmt = '%16.6f';
in.fid = fid;
mprint(tmp,in);

xibmean=mean(xbeta);
xmean=mean(x);
olsba=zeros(neq,2);
varba=zeros(neq,2);
tba=zeros(neq,2);
for i=1:m
    olsba(i,1)=xibmean*par((i-1)*nobeq+1:i*nobeq);
    olsba(i,2)=xmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    for j1=1:nobeq
        for j2=1:nobeq
            varba(i,1)=varba(i,1)+xibmean(j1)*xibmean(j2)*xpxi((i-1)*nobeq+j1,(i-1)*nobeq+j2);
        end
    end
    for j1=1:npar(i)
        for j2=1:npar(i)
            varba(i,2)=varba(i,2)+xmean(cumpar(i)+j1)*xmean(cumpar(i)+j2)*xpxi(m*nobeq+cumpar(i)+j1,m*nobeq+cumpar(i)+j2);
        end
    end
    if (varba(i,1)>0) tba(i,1)=olsba(i,1)/sqrt(varba(i,1)); end
    if (varba(i,2)>0) tba(i,2)=olsba(i,2)/sqrt(varba(i,2)); end
end
for j=1:nobeq
    huulp=0;
    for i=1:m
        huulp=huulp+par((i-1)*nobeq+j);
    end
    if (j==1) parhulp(j,1)=1-huulp;else parhulp(j,1)=-huulp; end
end
for i=neq
    olsba(i,1)=xibmean*parhulp(1:nobeq);
    olsba(i,2)=xmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    for j1=1:nobeq
        for j2=1:nobeq
            for j=1:m
                for p=1:m
                    varba(i,1)=varba(i,1)+xibmean(j1)*xibmean(j2)*xpxi((j-1)*nobeq+j1,(p-1)*nobeq+j2);
                end
            end
        end
    end
    for j1=1:npar(i)
        for j2=1:npar(i)
            varba(i,2)=varba(i,2)+xmean(cumpar(i)+j1)*xmean(cumpar(i)+j2)*xpxi(m*nobeq+cumpar(i)+j1,m*nobeq+cumpar(i)+j2);
        end
    end
    if (varba(i,1)>0) tba(i,1)=olsba(i,1)/sqrt(varba(i,1)); end
    if (varba(i,2)>0) tba(i,2)=olsba(i,2)/sqrt(varba(i,2)); end
end
[olsba(:,1) tba(:,1) olsba(:,2) tba(:,2)]

sumpar=zeros(nobeq,1);
for j=1:m
    j1=(j-1)*nobeq+1;
    j2=j*nobeq;
    sumpar=sumpar+par(j1:j2);
end
betafactor=zeros(nobs,neq);
for j=1:neq-1
    j1=(j-1)*nobeq+1;
    j2=j*nobeq;
    for i=1:nobs
        betafactor(i,j)=xbeta(i,:)*par((j-1)*nobeq+1:j*nobeq);
    end
end
j=neq;
for i=1:nobs
    betafactor(i,j)=xbeta(i,1)*(1-sumpar(1))-xbeta(i,2:end)*sumpar(2:nobeq);
end
sigmavec=zeros(0.5*m*(m+1),1);
tel=0;
for i=1:m
    for j=i:m
        tel=tel+1;
        sigmavec(tel)=sigma(i,j);
    end
end
ML_lesnonspat=-f2_lessarquotient([par;zeros(neq,1);sigmavec],nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,W)
% ***********************End estimation of LES without spatial effects ******
%
% Estimation with spatial effects and cost differences
%
global functel;
functel=0;
parold=[par;deltasar];
options = optimset('fminsearch');
options.MaxFunEvals=1000;
clear yhulp xhulp ymat itertot par;
iter=0;
converge=1.0;
criteria=0.001;
itermax=25; %
ptel=functel;
%
tic
while ( converge>criteria && iter<itermax)
iter=iter+1
par=fminsearch('f_lessarquotient',parold,options,nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,sigmai,W);
[functel-ptel f_lessarquotient(par,nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,sigmai,W)]
ptel=functel;
beta=par(1:nobeta);
alpha=par(nobeta+1:nobeta+nogam);
delta=par(nobeta+nogam+1:nobeta+nogam+neq);
betafactor=zeros(nobs,neq);
alphaf=zeros(nobs,neq);
alphabreuk=zeros(nobs,neq);
res=zeros(nobs,m);
hulp=zeros(nobs,1);

sumpar=zeros(nobeq,1);
for j=1:m
    sumpar=sumpar+beta((j-1)*nobeq+1:j*nobeq);
    betafactor(:,j)=xbeta*beta((j-1)*nobeq+1:j*nobeq);
end
betafactor(:,neq)=xbeta(:,1)*(1-sumpar(1))-xbeta(:,2:end)*sumpar(2:nobeq);

for j=1:neq
    alphaf(:,j)=x(:,1+cumpar(j):cumpar(j+1))*alpha(1+cumpar(j):cumpar(j+1));
    alphabreuk(:,j)=alphaf(:,j)./(wx(:,1+cumpar(j):cumpar(j+1))*alpha(1+cumpar(j):cumpar(j+1)));
end

for i=1:nobs
    for j=1:neq
        hulp(i)=hulp(i)+delta(j)*wy(i,j)*alphabreuk(i,j)+alphaf(i,j);
    end
    for j=1:m
        res(i,j)=y(i,j)-delta(j)*wy(i,j)*alphabreuk(i,j)-alphaf(i,j)-betafactor(i,j)*(xi(i)-hulp(i));
    end
end

for i=1:m
 for j=i:m
    sigma(i,j) = (res(:,i)'*res(:,j))/nobs;
    if j > i; sigma(j,i) = sigma(i,j); end
end
end
sigmai=inv(sigma);

parspit(:,iter)=par;
converge=sum(abs(par-parold));
parold=par;

save results; % save results in case program needs to be restarted later 
end %while loop
toc
iter
%
% Measures of R-squared
%
for i=1:nobs
    ymat(i,:)=y(i,:)-mean(y);
end
emat=zeros(nobs,neq);
for i=1:nobs
    for j=1:neq
        emat(i,j)=y(i,j)-delta(j)*wy(i,j)*alphabreuk(i,j)-alphaf(i,j)-betafactor(i,j)*(xi(i)-hulp(i));
    end
end
xt=emat'*emat;
xn=ymat'*ymat;
for i=1:neq
    R2(i)=1-xt(i,i)/xn(i,i);
end
R2sys2=R2*diag(xt)/sum(diag(xt));

xtel=0;
xnoem=0;
for i=1:nobs
    xtel=xtel+emat(i,1:m)*inv(sigma(1:m,1:m))*emat(i,1:m)';
    xnoem=xnoem+ymat(i,1:m)*inv(sigma(1:m,1:m))*ymat(i,1:m)';
end
R2sys1=1-xtel/xnoem;
% +++++++++++++++++++++
%
% Calculation var-cov matrix for T-values
%
% +++++++++++++++++++++
sigmavec=zeros(0.5*m*(m+1),1);
tel=0;
for i=1:m
    for j=i:m
        tel=tel+1;
        sigmavec(tel)=sigma(i,j);
    end
end
parm=[par;sigmavec];
tic
dhessn = hessian('f2_lessarquotient',parm,nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,W);
toc
xpxi = invpd(dhessn);
tvar = abs(diag(xpxi));
tstat = par./sqrt(tvar(1:nobeta+nogam+neq,1));
d1names=strvcat('delta1','delta2','delta3','delta4','delta5','delta6','delta7');
vnames=strvcat(vLES,v1beta,v2beta,v3beta,v4beta,v5beta,v6beta,v1names,v2names,v3names,v4names,v5names,v6names,v7names,d1names); % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7
%print results
fid=1;
fprintf(fid,'\n');
for i=1:neq
    fprintf(fid,'R-squared eq.   = %2d,%9.4f   \n',i,R2(i));
end
fprintf(fid,'R-squared system McElroy en Dhrymes  = %11.4f,%11.4f  \n',R2sys1,R2sys2);
% now print coefficient estimates, t-statistics and probabilities
tout = norm_prb(tstat); % find asymptotic z (normal) probabilities
tmp = [par tstat tout];  % matrix to be printed
% column labels for printing results
bstring = 'Coefficient'; tstring = 'Asymptot t-stat'; pstring = 'z-probability';
cnames = strvcat(bstring,tstring,pstring);
in.cnames = cnames;
in.rnames = vnames;
in.fmt = '%16.6f';
in.fid = fid;
mprint(tmp,in);
fprintf(fid,'Coefficients and T-values betas of equation that was eliminated in the estimation due to adding-up restrictions\n');
betaM=zeros(nobeq,1);
varM=zeros(nobeq,1);
for k=1:nobeq
    for j=1:m
        betaM(k)=+betaM(k)-par((j-1)*nobeq+k);
        for p=1:m
            varM(k)=varM(k)+xpxi((j-1)*nobeq+k,(p-1)*nobeq+k);
        end
    end
end
betaM(1)=betaM(1)+1;
[betaM betaM./sqrt(varM)]
xibmean=mean(xbeta);
xmean=mean(x);
wxmean=mean(wx);
wymean=mean(wy);
olsba=zeros(neq,3);
varba=zeros(neq,2);
olsce=zeros(neq,1);
varce=zeros(neq,1);
tba=zeros(neq,3);
for i=1:m
    olsba(i,1)=xibmean*par((i-1)*nobeq+1:i*nobeq);
    olsba(i,2)=xmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    olsba(i,3)=wxmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    olsce(i,1)=olsba(i,2)/olsba(i,3)*delta(i)*wymean(i)+olsba(i,2);
    for j1=1:nobeq
        for j2=1:nobeq
            varba(i,1)=varba(i,1)+xibmean(j1)*xibmean(j2)*xpxi((i-1)*nobeq+j1,(i-1)*nobeq+j2);
        end
    end
    for j1=1:npar(i)
        for j2=1:npar(i)
            varba(i,2)=varba(i,2)+xmean(cumpar(i)+j1)*xmean(cumpar(i)+j2)*xpxi(m*nobeq+cumpar(i)+j1,m*nobeq+cumpar(i)+j2);
        end
    end
    if (varba(i,1)>0) tba(i,1)=olsba(i,1)/sqrt(varba(i,1)); end
    if (varba(i,2)>0) tba(i,2)=olsba(i,2)/sqrt(varba(i,2)); end
end
for j=1:nobeq
    huulp=0;
    for i=1:m
        huulp=huulp+par((i-1)*nobeq+j);
    end
    if (j==1) parhulp(j,1)=1-huulp;else parhulp(j,1)=-huulp; end
end
for i=neq
    olsba(i,1)=xibmean*parhulp(1:nobeq);
    olsba(i,2)=xmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    olsba(i,3)=wxmean(1,cumpar(i)+1:cumpar(i+1))*par(m*nobeq+cumpar(i)+1:m*nobeq+cumpar(i+1));
    olsce(i,1)=olsba(i,2)/olsba(i,3)*delta(i)*wymean(i)+olsba(i,2);
    for j1=1:nobeq
        for j2=1:nobeq
            for j=1:m
                for p=1:m
                    varba(i,1)=varba(i,1)+xibmean(j1)*xibmean(j2)*xpxi((j-1)*nobeq+j1,(p-1)*nobeq+j2);
                end
            end
        end
    end
    for j1=1:npar(i)
        for j2=1:npar(i)
            varba(i,2)=varba(i,2)+xmean(cumpar(i)+j1)*xmean(cumpar(i)+j2)*xpxi(m*nobeq+cumpar(i)+j1,m*nobeq+cumpar(i)+j2);
        end
    end
    if (varba(i,1)>0) tba(i,1)=olsba(i,1)/sqrt(varba(i,1)); end
    if (varba(i,2)>0) tba(i,2)=olsba(i,2)/sqrt(varba(i,2)); end
end
df1=zeros(neq,length(par)-m*nobeq);
for i=1:neq
    df1(i,cumpar(i)+1:cumpar(i+1))=(xmean(1,cumpar(i)+1:cumpar(i+1))*olsba(i,3)-wxmean(1,cumpar(i)+1:cumpar(i+1))*olsba(i,2))/...
        (olsba(i,3)^2)*delta(i)*wymean(1,i)+olsba(i,2);
    df1(i,nogam+i)=olsba(i,2)/olsba(i,3)*wymean(1,i);
    varce(i,1)=df1(i,:)*xpxi(m*nobeq+1:m*nobeq+nogam+neq,m*nobeq+1:m*nobeq+nogam+neq)*df1(i,:)';
    if (varce(i,1)>0) tba(i,3)=olsce(i,1)/sqrt(varce(i,1)); end
end
tba=real(tba);
[olsba(:,1) tba(:,1) olsba(:,2) tba(:,2) olsce(:,1) tba(:,3)]
loglik=-f2_lessarquotient(parm,nobs,neq,m,nogam,nobeta,nobeq,cumpar,y,wy,x,wx,xi,xbeta,W)
% Test hypothesis deltas LES = delta single equations
R=eye(neq);
Rs=R*delta-deltasar;
Wald=Rs'*inv(R*xpxi(totpar+1:totpar+neq,totpar+1:totpar+neq)*R')*Rs/neq
1-chis_cdf(real(Wald)*neq,neq) % probability greater than 0.05 points to insignificance
1-fdis_cdf(real(Wald),neq,nobs-totpar-neq)% probability greater than 0.05 points to insignificance
R=zeros(m,m*nobeq);
for i=1:m
    for j=1:nobeq
        R(i,(i-1)*nobeq+j)=xibmean(j);
    end
end
% Test hypothesis betas LES = betas LES without interaction effects
Rs=R*par(1:m*nobeq)-[0.0653;0.0473;0.3014;0.1619;0.1084;0.2247]; % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7, 
% I filled in the deltas LES without interaction effects from a previous run
% of this program
Wald=Rs'*inv(R*xpxi(1:m*nobeq,1:m*nobeq)*R')*Rs/m
1-chis_cdf(real(Wald)*m, m) % probability greater than 0.05 points to insignificance
1-fdis_cdf(real(Wald),m,nobs-totpar-neq) % probability greater than 0.05 points to insignificance
R=zeros(neq,totpar-m*nobeq);
for i=1:neq
    for j=1:npar(i)
        R(i,cumpar(i)+j)=xmean(1,cumpar(i)+j);
    end
end
% Test hypothesis alphas LES = alphas LES without interaction effects
Rs=R*par(m*nobeq+1:totpar)-[-0.1438;0.0592;0.2179;0.1639;0.0669;0.1279;0.0699]; % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7, 
% I filled in the alphas LES without interaction effects from a previous run
% of this program
Wald=Rs'*inv(R*xpxi(m*nobeq+1:totpar,m*nobeq+1:totpar)*R')*Rs/neq
1-chis_cdf(real(Wald)*neq,neq) % probability greater than 0.05 points to insignificance
1-fdis_cdf(real(Wald),neq,nobs-totpar-neq)% probability greater than 0.05 points to insignificance
% Test hypothesis committed expenditures LES = committed expenditures LES without interaction effects
Rs=olsce(:,1)-[-0.1438;0.0592;0.2179;0.1639;0.0669;0.1279;0.0699]; % NOTE: NEEDS TO BE ADJUSTED IF neq IS NOT EQUAL TO 7, 
% I filled in the betas LES without interaction effects from a previous run
% of this program
Wald=Rs'*inv(df1*xpxi(m*nobeq+1:m*nobeq+nogam+neq,m*nobeq+1:m*nobeq+nogam+neq)*df1')*Rs/neq
1-chis_cdf(real(Wald)*neq,neq) % probability greater than 0.05 points to insignificance
1-fdis_cdf(real(Wald),neq,nobs-totpar-neq)% probability greater than 0.05 points to insignificance