function [simout,cmout]=runfomc(bindef,data,piPrior,PPrior,outfile)
%function [simout,cmout]=runfomc(bindef,data,piPrior,PPrior,outfile)
%This function uses BACC to estimate FOMC. Used to produce output for 
%"Inter-State Dynamics of Invention Activities, 1930--2000" joint with
%Catherine Co and Myeong-Su Yun.
%Inputs:
%       bindef: nbin x 2 vector with first column equal to lower bin
%       definition and second column equal to upper bin definition
%       data: T x N matrix of raw data where T equals number of time series
%       observations and N is number of cross sectional units
%       piPrior: (1 x nbin) prior matrix for initial distirbution
%       PPrior: (nbin x nbin) prior matrix for transition matrix
%       outfile: name of file to print results to
%Outputs:
%       simout: raw simulation output structure;
%       cmout: raw concentration measure output structure;
% Author: John Landon-Lane
%         lane@econ.rutgers.edu
%         75 Hamilton Street
%         New Brunswick, NJ 08901
%         USA
% Date of last revision: Oct 30, 2005
% Notes: each column of data matrix is a year. For the 10-year transition
% case the years are 1930, 1940, 1950, 1960, 1970, 1980, 1990, and 2000. It
% is assumed that the user has created a data matrix with 8 rows with each
% row equally a year before inputting into program.

fid=fopen(outfile,'w');
[n,m]=size(PPrior);
nsim=10000;
nbin=size(bindef,1);



% set up data matrix S to input into BACC routine

S=Smat(bindef,data);


%set up model instance (BACC)
mi=minst('nsfomfs','pi','P',piPrior,PPrior,S); %BACC command
sim=extract(mi); % BACC command

%draw from prior
sim=priorsim(mi,nsim,1); %BACC command
%draw from posterior
sim=postsim(mi,nsim,1); %BACC command


%calculate prior moments

[piprmn,piprse]=dirmoment(piPrior); % prior moments for initial distribution
[Pprmn,Pprse]=dirmoment(PPrior); % prior moments for transition matrix
invpr=invget(sim.PPrior); % invariant disitrbution for prior draws
invprmn=mean(invpr,3); % mean of invariant distribution implied by prior
invprse=std(invpr,1,3); % standard error of invariant disitribution implied by prior

%print results to file
fprintf(fid,'Prior Moments\n');
fprintf(fid,'-------------\n\n');

fprintf(fid,'Initial Distribution\n');
for i=1:n
    fprintf(fid,' %4.3f   ',piprmn(i));
end
fprintf(fid,'\n');
for i=1:n
    fprintf(fid,'(%4.3f)  ',piprse(i));
end
fprintf(fid,'\n\n');

fprintf(fid,'Transition Matrix\n');
for j=1:n
    for i=1:n
        fprintf(fid,' %4.3f   ',Pprmn(j,i));
    end
    fprintf(fid,'\n');
    for i=1:n
        fprintf(fid,'(%4.3f)  ',Pprse(j,i));
    end
    fprintf(fid,'\n');
end
fprintf(fid,'\n\n');
        
fprintf(fid,'Invariant Distribution\n');
for i=1:n
    fprintf(fid,' %4.3f   ',invprmn(i));
end
fprintf(fid,'\n');
for i=1:n
    fprintf(fid,'(%4.3f)  ',invprse(i));
end
fprintf(fid,'\n\n');

%calculate posterior moments
pimn=mean(sim.pi,3); %mean posterior for initial distribution
pise=std(sim.pi,1,3); % posterior SE for initial distribution
Pmn=mean(sim.P,3); % posterior mean of transition matrix
Pse=std(sim.P,1,3); % posterior SE of transition matrix

%calulate invariant distribution
inv=invget(sim.P); 
invmn=mean(inv,3); % posterior mean of invariant distribution
invse=std(inv,1,3); % posterior SE of invariant distirbution

% print output to file
fprintf(fid,'Posterior Moments\n');
fprintf(fid,'-------------\n\n');

fprintf(fid,'Initial Distribution\n');
for i=1:n
    fprintf(fid,' %4.3f   ',pimn(i));
end
fprintf(fid,'\n');
for i=1:n
    fprintf(fid,'(%4.3f)  ',pise(i));
end
fprintf(fid,'\n\n');

fprintf(fid,'Transition Matrix\n');
for j=1:n
    for i=1:n
        fprintf(fid,' %4.3f   ',Pmn(j,i));
    end
    fprintf(fid,'\n');
    for i=1:n
        fprintf(fid,'(%4.3f)  ',Pse(j,i));
    end
    fprintf(fid,'\n');
end
fprintf(fid,'\n\n');
        
fprintf(fid,'Invariant Distribution\n');
for i=1:n
    fprintf(fid,' %4.3f   ',invmn(i));
end
fprintf(fid,'\n');
for i=1:n
    fprintf(fid,'(%4.3f)  ',invse(i));
end
fprintf(fid,'\n\n');



% calculate concentration measure for intiial and invariant distributions

cmmean0=conmean(sim.pi); % standard deviation of initial distribution
cmmeani=conmean(inv); % standard deviation of invariant distribution

cmout.sd0=cmmean0; % output 
cmout.sdi=cmmeani; % output
mcmm0=mean(cmmean0);
scmm0=std(cmmean0,1);

mcmmi=mean(cmmeani);
scmmi=std(cmmeani,1);

count=0;
for i=1:nsim
    if(cmmeani(i)<cmmean0(i))
        count=count+1;
    end
end
%print results to file
fprintf(fid,'\n\n');
fprintf(fid,'Std Dev Initial Distirbution\n');
fprintf(fid,' %6.4f \n',mcmm0);
fprintf(fid,'(%6.4f)\n',scmm0);

fprintf(fid,'\n\n');
fprintf(fid,'Std Dev limiting dist\n');
fprintf(fid,' %6.4f \n',mcmmi);
fprintf(fid,'(%6.4f)\n',scmmi);


fprintf(fid,'\n\n');
fprintf(fid,'Posterior Probability that SD(inv)<SD(init)\n');
fprintf(fid,'%4.3f',count/nsim);

%Calculate hpd of difference between pi0 and inv

diff=cmmeani-cmmean0;

mdiff=mean(diff);
sddiff=std(diff,1);

fprintf(fid,'\n\n');
fprintf(fid,'Difference in Std Dev limiting dist\n');
fprintf(fid,' %6.4f \n',mdiff);
fprintf(fid,'(%6.4f)\n',sddiff);


p=[0.01:0.01:0.7]';
np=size(p,1);
for i=1:np
    hpddiff(i,:)=hpd(diff,p(i));
    hpdinit(i,:)=hpd(cmmean0,p(i));
    hpdlim(i,:)=hpd(cmmeani,p(i));
end


fprintf(fid,'\n\n');
fprintf(fid,'HPD: Std Dev of distribution Initial\n');
for i=1:np
    fprintf(fid,'%2.0f percent HPD: [%4.3f, %4.3f]\n',100*(1-p(i)),hpdinit(i,:));
end 

fprintf(fid,'\n\n');
fprintf(fid,'HPD: Std Dev of distribution: Limiting\n');
for i=1:np
    fprintf(fid,'%2.0f percent HPD: [%4.3f, %4.3f]\n',100*(1-p(i)),hpdlim(i,:));
end 

fprintf(fid,'\n\n');
fprintf(fid,'HPD: Std Dev of distribution: Difference\n');
for i=1:np
    fprintf(fid,'%2.0f percent HPD: [%4.3f, %4.3f]\n',100*(1-p(i)),hpddiff(i,:));
end    

simout=sim; %output file

% mobility measures

mob=mobcalc(sim.P);
mobmn=mean(mob,3);
mobsd=std(mob,1,3);
M=n;
% print out mobility indices
fprintf(fid,'\n\n');
fprintf(fid,'Mobility Indices\n');
fprintf(fid,'----------------\n');
fprintf(fid,'Shorrocks Overall Measure: %4.3f\n',mobmn(1,1));
fprintf(fid,'                          (%4.3f)\n\n',mobsd(1,1));
fprintf(fid,'Upward Overall Measure: %4.3f\n',mobmn(1,2));
fprintf(fid,'                       (%4.3f)\n\n',mobsd(1,2));
fprintf(fid,'Downward Overall Measure: %4.3f\n',mobmn(1,3));
fprintf(fid,'                         (%4.3f)\n\n',mobsd(1,3));
for i=1:M
    fprintf(fid,'Prais measure class %2.0f: %4.3f\n',i,mobmn(1,3+i));
    fprintf(fid,'                       (%4.3f)\n\n',mobsd(1,3+i));
end
for i=1:M-1
    fprintf(fid,'Prob of moving up (class %2.0f): %4.3f\n',i,mobmn(1,3+M+i));
    fprintf(fid,'                             (%4.3f)\n\n',mobsd(1,3+M+i));
end    
for i=2:M
    fprintf(fid,'Prob of moving down (class %2.0f): %4.3f\n',i,mobmn(1,3+M+(M-1)+(i-1)));
    fprintf(fid,'                               (%4.3f)\n\n',mobsd(1,3+M+(M-1)+(i-1)));
end
fprintf(fid,'Modulus of second eigenvalue: %4.3f\n',mobmn(1,3+M+2*(M-1)+1));
fprintf(fid,'                             (%4.3f)\n',mobsd(1,3+M+2*(M-1)+1));
fprintf(fid,'Half Life: %4.3f\n',mobmn(1,3+M+2*(M-1)+2));
fprintf(fid,'          (%4.3f)\n',mobsd(1,3+M+2*(M-1)+2));


fclose(fid);

%--------------------------------------------------------------------------
function cmout=conmean(xin)
%THis function calcualtes the std dev of a discrete distirbution
[n1,n2,n3]=size(xin);
bin=[1:n2];
for i=1:n3
    %calcualte mean of iteration i
    mn=sum(bin.*xin(1,:,i));
    %calculate std of iteration i
    b2=(bin-mn).^2;
    var=sum(b2.*xin(1,:,i));
    cmout(i)=sqrt(var);
end


%--------------------------------------------------------------------------
function mout=mobcalc(Pin)
%Calculates set of mobility measures

% mout is 1xnxnsim array
% mobility measures calculates are
%1. overall shorrocks
%2. upward SHorrocks
%3. downward Shorrocks
%4. prob of moving (M of them)
%5. upward Prais (M-1) of them
%6. downward Prais (M-1) of them
%7. modulus of second eigenvalue

[m1,m2,nsim]=size(Pin);

if(m1~=m2)
    error('Input array not square transition matrix');
end
nout=3+m1+2*(m1-1)+2;
mout=zeros(1,nout,nsim);
M=m1;
for i=1:nsim
    P=Pin(:,:,i);
    mout(1,1,i)=(M-trace(P))/(M-1);
    for j=1:M
        mout(1,3+j,i)=(1-P(j,j));
    end
    for j=1:M-1
        mout(1,3+M+j,i)=sum(P(j,[j+1:M]));
    end
    for j=2:M
        mout(1,3+M+(M-1)+(j-1),i)=sum(P(j,[1:j-1]));
    end
    mout(1,2,i)=sum(mout(1,[3+M+1:3+M+M-1],i))/(M-1);
    mout(1,3,i)=sum(mout(1,[3+M+(M-1)+1:3+M+(M-1)+M-1],i))/(M-1);
    e=eig(P);
    e=sort(e);
    mout(1,nout-1,i)=abs(e(M-1));
    mout(1,nout,i)=-log(2)/log(abs(e(M-1)));
end
%--------------------------------------------------------------------------
function Sout=Smat(bdef,datin);
%Calculates S mat from data using bin definitions --used as input to BACC
%routines
%Sout=Smat(bdef,datin)

[T,N]=size(datin);

[nb,mb]=size(bdef);


% calculation of sout matrix
for i=1:T
    for k=1:N
        for j=1:nb
            if(datin(i,k) > bdef(j,1))
                if(datin(i,k) <= bdef(j,2))
                    Sout(i,k)=j;
                end
            end
        end
    end
end
%--------------------------------------------------------------------------
function inv=invget(P)
%Computes invariant distribution for transition matrices given in P

[n1,n2,nsim]=size(P);
if(n1~=n2)
    error('Input to INVGET is not a transition matrix');
end

for i=1:nsim
    pm=P(:,:,i);
   [w,d]=eig(pm.');
   w=conj(w);
   w=real(w/sum(w));
   inv(1,:,i)=w(:,1)';
end
%--------------------------------------------------------------------------
function [mean,std]=dirmoment(a)
%This function takes in parameters defining a Dirichlet RV and outputs the
%theoretical mean and standard deviation

[n,m]=size(a);

sa=sum(a,2);

mean=a./kron(sa,ones(1,m));

ee=mean.*(1-mean);
var=ee./kron(1+sa,ones(1,m));
std=var.^(0.5);
%--------------------------------------------------------------------------
function [hpdout]=hpd(xin,p,nsteps)
% [hi,lo]=hpd(xin,p)
%This routine returns the (1-p)% highest posterior density region for a
%vector input xin.
%
%Inputs:
%       xin: vector of observations
%       p: probability
%Outputs:
%       hi: upper bound
%       lo: lower bound of hpd

[f,x]=ecdf(xin);

lo=min(xin);
hi=max(xin);

[n,m]=size(f);

for i=1:n
    if(f(i)<=p/2)
        lind=i;
    end
    if(f(i)<=(1-p/2))
        uind=i;
    end
end
lo=x(lind);
hi=x(uind);

hpdout=[lo hi];

