%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finite Mixture Regression Model (FMRM)
% based on KPT section 15.6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;  % start stop watch

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

load c:\klaus\AAEC6564\mlab\worksp\mod10_FMRM_data_v2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=length(y);
G=3; %number of classes
k=size(X,2); %we use the same X for all regimes, so only 1 k is needed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=20000; % burn-ins
r2=10000; % keepers
R=r1+r2;

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(n-k);

% elements for beta (same prior for all betas)
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta
betadraw=bols; %use OLS values as starting draws for Gibbs Sampler
%same for all betas

% elements for sig2 (same prior for all sigmas)
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
%same for all sigma2's
v0=1/2;
tau0=1/2; % diffuse prior shape and scale

%prior for cp
acp0=[2;2;2]; %used by KPT - don't need to sum to 1

%draw starting vector of cp
cpdraw=mean(dirrnd(acp0,100),2); %should be approx. 1/3 each

%draw starting vector of indicators
Zdraw=mnrnd(1,cpdraw,n); %n by 3

% break
% 'ok'

[beta1mat,beta2mat,beta3mat,sig11mat,sig22mat,sig33mat,cpmat,Zmat]=...
    gs_fmrm(X,y,k,n,G,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw,...
    acp0,cpdraw,Zdraw);


'GS done'
% put all draws together & run diagnostics
allmat=[beta1mat;beta2mat;beta3mat;sig11mat;sig22mat;sig33mat;cpmat];
kdiag=klausdiagnostics_greater0(allmat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Capture unordered results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Open log file 
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod10_fmrm_v2.txt','w');
if fid==-1;
   warning('File could not be opened');
else;
   disp('File opened successfully');
end;


fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');

fprintf(fid,'true beta1\t%6.3f\n',beta1true'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'true beta2\t%6.3f\n',beta2true'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'true beta3\t%6.3f\n',beta3true'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

sigAlltrue=[sig1true;sig2true;sig3true];
fprintf(fid,'true sig2s \t%6.3f\n',sigAlltrue'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

fprintf(fid,'true regime probs \t%6.3f\n',cp'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

% beta stuff
out=kdiag(1:k,:)';
fprintf(fid,'Output table for beta1 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(k+1:2*k,:)';
fprintf(fid,'Output table for beta2 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(2*k+1:3*k,:)';
fprintf(fid,'Output table for beta3 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(3*k+1:3*k+3,:)';
fprintf(fid,'Output table for sig2s \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(3*k+4:end,:)';
fprintf(fid,'Output table for cp \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


save c:\klaus\AAEC6564\mlab\worksp\mod10_fmrm_v2 beta1mat beta2mat beta3mat ...
    sig11mat sig22mat sig33mat cpmat Zmat;


finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);

st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;


