% CLUSTER ANALYSIS RESULTS
% March 26, 2010 12:30pm

% LOAD DATA %
%%%%%%%%%%%%%

% Data generated by mmp_cluster_data9.dta - RE-RUN
clear
cd('/Users/fgarip/Documents/Clustering Paper 07/Analysis/Matlab Analysis');
importfile('/Users/fgarip/Documents/Data/MMP Data/Clustering/rdata1.raw');

% cd('M:\Documents\C Manuscripts\Clustering Paper 07\Analysis\Matlab Analysis');
% importfile('M:\Documents\Data\MMP Data\Clustering\rdata1.raw');

seed = 931316785;
rand('seed',seed);
randn('seed',seed);

head = (relhead==1);
met = (metro<=2);

id = [commun hhnum];
id = id(find(head==1),:);

usemp = usemp/100000;
lnfdi = log(fdi);
lw = lw/100000;

    % Continuous versions of variables used in clustering

datac = [head sex educ ocagri ocmanuf ocserv room lnland business mxmig nhhleg ...
         nhhnleg pcmig agri self ltmin lnpop];

datanamesc = {'head','sex','educ','ocagri','ocmanuf','ocserv',...
            'room','lnland','business','mxmig','nhhleg','nhhnleg','pcmig',...
            'agri','self','ltmin','lnpop'};
        
    % Variables not used in clustering

datao = [rain km prisch secsch bank road min2 dagri ejido manu dmanu lnmanp ...
         lfp dlfp inf mxint mxwage uswage wratio lnexports lnimports lntrade ...
         usunemp usemp usempgr lnfdi fdigr lnexports lnimports lntrade ...
         tradegr lw app visa lnins lnbpe ];

datanameso = {'rain','km','prisch','secsch','bank','road','min2','dagri', ...
              'ejido','manu','dmanu','lnmanp','lfp','dlfp','inf','mxint', ...
              'mxwage','uswage','wratio','exports','imports','trade', ...
              'usunemp','usemp','usempgr','lnfdi','fdigr','lnexports', ... 
              'lnimports','lntrade','tradegr','lw','app','visa','lnins','lnbpe'};
   
    % Variables used in clustering - ALL BINARY, N=17

data = [head sex sec sec2 ocagri ocmanuf room1 land1 business1 ...
        mxmig hhleg hhnleg pcmig1 agri1 self1 ltmin1 met];
   
datanames = {'head','sex','sec','sec2','ocagri','ocmanuf',...
            'room1','land1','business1','mxmig','hhleg','hhnleg','pcmig1',...
            'agri1','self1','ltmin1','met'};

datacur = data;

% K-MEANS CLUSTERING %
%%%%%%%%%%%%%%%%%%%%%%
   
       % IMPORTANT NOTE: We need to run many replicates for kmeans. Matlab
       % automatically selects the solution that minimizes the total sum of
       % distances (see the demo of cluster analysis). Otherwise, we may
       % get stuck in local minima, and miss the global solution.
        
% [cl2,cm2] = kmeans(datacur,2, 'distance','cityblock','replicates',50,'display','final'); 
% [s2,h2] = silhouette(datacur, cl2,'cityblock');
% saveas(h2,'sil2.tif');
 
%[cl3,cm3] = kmeans(datacur,3,'distance','cityblock','replicates',100,'display','final');
%[cl4,cm4] = kmeans(datacur,4,'distance','cityblock','replicates',1000,'display','final'); 

% Select the number of clusters

%nocl = 4;
%clcur = cl4;
%cmcur = cm4;

   % Save the clusters for hh heads only for regression analysis later on

    %clcur_head = clcur(head==1);

%save CLUSTERS data datanames cl* cm* nocl id


load CLUSTERS

% Re-define data and datanames vectors to preserve the order of mxmig (it
% is different in saved CLUSTERS results)

data = [head sex sec sec2 ocagri ocmanuf room1 land1 business1 ...
        mxmig hhleg hhnleg pcmig1 agri1 self1 ltmin1 met];
   
datanames = {'head','sex','sec','sec2','ocagri','ocmanuf',...
            'room1','land1','business1','mxmig','hhleg','hhnleg','pcmig1',...
            'agri1','self1','ltmin1','met'};


    % Re-order clusters such that cluster 1: Income maximizers
    % 2: Rural risk diversifiers, 3: Network migrants, 
    % 4: New urban stream

    cltemp = clcur;
    cltemp(clcur==1) = 1;
    cltemp(clcur==3) = 2;
    cltemp(clcur==4) = 3;
    cltemp(clcur==2) = 4;
    
    clcur = cltemp;
    clear cltemp;
    

    % Transfer data to R for heatmaps
    
    % Outcomes of interest
    % 61 ttrip obs missing - set to zero
    
    mpertrip = texp./ttrip;
    mpertrip(isnan(mpertrip)==1)=0;
    
    reg = [mpertrip repmig];
       
     
        % Compute distance of each observation to cluster centroids. This
        % is necessary for some of the heatmap analysis in R - i.e., when
        % we need to select a small number of observations, we select those
        % that are closest to centroids.
     
     dist_centroid = zeros(size(clcur,1), max(clcur));

     tic
     for i=1:size(dist_centroid,1)
         for j=1:size(dist_centroid,2)
             dist_centroid(i,j) = sum(abs(data(i,:)-cmcur(j,:)));
         end
     end
     toc
     
     dist = zeros(size(clcur,1),1);
     for j=1:max(clcur) 
        dist(clcur==j) = dist_centroid(clcur==j,j);
     end

    cd('/Users/fgarip/Documents/Clustering Paper 07/Analysis/');
    %cd('M:\Documents\C Manuscripts\Clustering Paper 07\Analysis\');

    dlmwrite('cluster.txt', clcur,',');
    dlmwrite('reg.txt',reg,',');
    dlmwrite('dist.txt',dist,',');
    
    % Write data to a text file for Mexico Maps
    % Create a matrix for each year (1970,1980,1990)[including the
    % following 5 years] that records the % of individuals in each 
    % of the 4 clusters in 21 states (21x4)

    uniq_state = unique(statenum);
    state_names = {'Aquascalientes'; 'Baja Norte'; 'Colima'; 'Chihuahua';...
        'Durango';'Guanajuato';'Guerrero';'Hidalgo';'Jalisco';'Mexico';'Michoacan';...
        'Morelos';'Nayarit';'Nuevo Leon'; 'Oaxaca'; 'Puebla';'San Luis Potosi';...
        'Sinaloa';'Tlaxcala';'Veracruz';'Zacatecas'};
    yrs = [1970;1980;1990];
    state_clust = zeros(size(uniq_state,1),max(clcur),size(yrs,1));
    state_clust_allyrs = zeros(size(uniq_state,1),max(clcur));
    
   
    for k = 1:size(yrs,1)
        for j=1:max(clcur)
            for i=1:size(uniq_state,1)
                state_clust(i,j,k) = sum(statenum==uniq_state(i) & clcur==j & year>=yrs(k) & year<=yrs(k)+5)*100/sum(statenum==uniq_state(i) & year>=yrs(k) & year<=yrs(k)+5);
            end
        end
    end
    
        % Combine all clusters observed only in the time period they were
        % dominant at. Cluster 1 in 1970-75, 2 in 1980-85, 3 and 4 both in 
        % 1990-95.
  
   state_clust_comb = [state_clust(:,1,1) state_clust(:,2,2) state_clust(:,3,3) state_clust(:,4,3)];
   
   
        % Compute cluster compositions in each state over the entire time
        % horizon
    
    for j=1:max(clcur)
        for i=1:size(uniq_state,1)
            state_clust_allyrs(i,j) = sum(statenum==uniq_state(i) & clcur==j)*100/sum(statenum==uniq_state(i));
        end
    end
    
        % Categorize the data: 1: 0-30% 2:31-50% 3: 51-100%
    
    state_clust(state_clust<=30)=1;
    state_clust(state_clust>30 & state_clust<=50)=2;
    state_clust(state_clust>50)=3;

    state_clust_7075 = state_clust(:,:,1);
    state_clust_8085 = state_clust(:,:,2);
    state_clust_9095 = state_clust(:,:,3);
    
    dlmwrite('state_nums.txt',uniq_state,',');
    dlmwrite('state_clust_comb.txt', state_clust_comb,',');
    dlmwrite('state_clust_allyrs.txt', state_clust_allyrs,',');
    dlmwrite('state_clust_7075.txt', state_clust_7075,',');
    dlmwrite('state_clust_8085.txt', state_clust_8085,',');
    dlmwrite('state_clust_9095.txt', state_clust_9095,',');
    
    % Record a slightly different statistic - % of income-maximizers from
    % each state in the entire time interval. This statistic may be biased since
    % different states are observed in each year. 
    
    clust_across_state = zeros(size(uniq_state,1),max(clcur));
    
    for j=1:max(clcur)
        for i=1:size(uniq_state,1)
            clust_across_state(i,j) = sum(statenum==uniq_state(i) & clcur==j)*100/sum(clcur==j);
        end
    end
    
    dlmwrite('clust_across_state.txt', clust_across_state,',');

% (A) Trends in Distribution of Migrants Across Clusters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
col = {'b-','r-','g-','y-','m-','k-'};

yr = [1970:1:2000]';
yclust = zeros(size(yr,1),nocl);

for i=1:size(yr)
    for j=1:nocl
        yclust(i,j) = sum(year(find(clcur==j))==1969+i)./sum(year==1969+i);
        yclustnum(i,j) = sum(year(find(clcur==j))==1969+i);
    end
end

    % Proportion of migrants %
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    
for j=1:nocl
    plot(yr,yclust(:,j),col{j})
    hold on
end
hold off

    % Number of migrants %
    %%%%%%%%%%%%%%%%%%%%%%

% for j=1:nocl
%     plot(yr,yclustnum(:,j),col{j})
%     hold on
% end
% hold off


    % Transfer data to R for plots
    
    yrclust = [yr yclust.*100];
    dlmwrite('yrclust.txt', yrclust,',');
    
    % Save annual trends keeping one observation per year
    % Transer to R for plots
    
    trends = [uswage inf visa lntrade year];
    trends = trends(year>=1970 & year<=2000,:);
    
    trends2 = [mxwage wratio mxint year];
    trends2 = trends2(year>=1970 & year<=2000,:);

    year2  = year(year>=1970 &  year<=2000);
    [srt, index] = unique(year2,'rows','first');
    trends  = trends(index,:);
    trends2 = trends2(index,:);
  
    dlmwrite('trends.txt', trends,',');
    dlmwrite('trends2.txt', trends2,',');
 
 

% (B) Cluster Centroid Characteristics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Means of variables used in clustering
% (binary and conts versions), and other variables

for j = 1:nocl
    means(:,j)      = mean(data(clcur==j,:))';
    meansc(:,j)     = mean(datac(clcur==j,:))';
    meanso(:,j)     = mean(datao(clcur==j,:))';
end


tabulate(clcur)

for i = 1:size(data,2)
    display(datanames(i));
    display('......................................');
    means(i,:)
end

for i = 1:size(datac,2)
    display(datanamesc(i));
    display('......................................');
    meansc(i,:)
end

for i = 1:size(datao,2)
    display(datanameso(i));
    display('......................................');
    meanso(i,:)
end


    % (B.1) T-test results %
    %%%%%%%%%%%%%%%%%%%%%%%%
    
    % Compare each cluster to the rest of the data
    % Ttest2 function assigns 1 if the two samples have different means
    % at 5% level, 0 o/w

test  = zeros(size(data,2),nocl);
testc = zeros(size(datac,2),nocl);
testo = zeros(size(datao,2),nocl);

for i = 1:size(data,2) 
    for j = 1:nocl
        test(i,j)      = ttest2(data(clcur==j,i), data(clcur~=j,i));
        testc(i,j)     = ttest2(datac(clcur==j,i), datac(clcur~=j,i));
        testo(i,j)     = ttest2(datao(clcur==j,i), datao(clcur~=j,i));
    end
end

    % Additional ttests for differentiating variables (shown in boldface on
    % table) - we compare the two closest clusters
    
    test_index = zeros(size(data,2),2);
    test_index = [1 4; 1 4; 4 3; 4 3; 1 2; 4 1; 2 4; 2 3; 2 4; 1 4; 3 2; 3 2; 3 2; 1 2; 2 1; 2 1; 4 3]
    test_add = zeros(size(data,2),1);
    
    for i = 1:size(data,2)
        test_add(i) = ttest2(data(clcur==test_index(i,1),i), data(clcur==test_index(i,2),i));
    end


    % Transfer numbers for the Excel tables

    tabulate(clcur)

    datanames'
    means
    test

    datanamesc'
    meansc
    testc

    datanameso'
    meanso 
    testo

    test_index
    test_add'
 
    

% (C) CHANGING REGRESSION COEFFICIENTS OVER TIME %
% MIGRANT-NONMIGRANT COMPARISON %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % Run separate figures by clusters, and show that the 
    % coefficients remain more stable over time.
    % That is, the change over time, is created by the different
    % types of clusters dominant at the time.
    
    % We use the life history data (only hh heads) - o/w non-migrants
    % are observed only on the survey year in the pers data.

cd('/Users/fgarip/Documents/Clustering Paper 07/Analysis/Matlab Analysis');
importfile('/Users/fgarip/Documents/Data/MMP Data/Clustering/fulldata.raw');

room1 = (troom>0);
land1 = (vland>0);
bus1  = (business>0);
sibusmig = (nosibusmig>0);
famusmig = (mousmiglag==1 | fausmiglag==1);

    
    % Load cluster membership information on hh heads, and match it with
    % the life history data
      
load CLUSTERS;

clust = zeros(size(mig,1),1);
for i = 1:size(id, 1)
    i
    ind = find(commun == id(i,1) & hhnum == id(i,2));
    clust(ind) = clcur_head(i);
end
 
    % There are a few mis-matches: 183 individuals have no cluster
    % membership, but seem to have migrated. 26 individuals seem in cluster
    % 2, without having migrated. 65 individuals seem in cluster 3,
    % without having migrated according to life history data. Assume that
    % pers data is more accurate. Set the former case to non-migrants, and
    % the latter to to migrants.

tabulate(mig(clust==0))
tabulate(mig(clust==1))
tabulate(mig(clust==2))
tabulate(mig(clust==3))
tabulate(mig(clust==4))

mig(clust==0) = 0;
mig(clust==2) = 1;
mig(clust==3) = 1;

datafull = [mig age sex sec sec2 mxmig ocagri ocmanuf room1 land1 bus1...
            famusmig prev agrim self ltmin met];
 

bg_year  = 1970;
en_year  = 2000;

data_int  = datafull(year>=bg_year & year<=en_year,:);
year_int  = year(year>=bg_year & year<=en_year);
clust_int = clust(year>=bg_year & year<=en_year);
id_int    = [commun hhnum persnum];
id_int    = id_int(year>=bg_year & year<=en_year,:);

y = data_int(:,1); 
x = data_int(:,2:end);
xnames = {'Age','Sex','Some Sec Education','Secondary Education',...
          'Mexican Migration','Agri. Occupation','Manuf. Occupation',...
          'Property Ownership','Land Ownership','Business Ownership',...
          'U.S. Mig. in Family','Mig. Prevalence',...
          '% in Agriculture','% Self Employed',...
          '% < Min Wage','Metropolitan'};

step = 1;
n = range(year_int)/step + 1;

coef = zeros(size(x,2)+1,n);

for t=1:n
     bg = bg_year + (t-1)*step;
     en = bg_year + t*step;    
     ind1 = find(year_int>=bg & year_int<en);
     
      % Keep non-migrants in the first period of the range

     [srt_id, ind2] = unique(id_int(ind1,:), 'rows', 'first');
     ind = ind1(ind2);
     
     coef(:,t)    = glmfit(x(ind,:),y(ind),'binomial','link','logit');
end

time = [bg_year:step:en_year];


for i = 1:size(coef,1)-1 
    figure(i);
    plot(time, coef(i+1,:),'k-');
    xlabel('{\bf Year}');
    ylabel(xnames{i});
    %close(i);
end


        % In an earlier analysis, I tried to show the stability of
        % coefficients for each migrant cluster over time. The results
        % suggested otherwise. We cannot expect coefficients similar within
        % a cluster over time as the proportion of migrants in that cluster
        % is changing over time as well. Since a small number of migrants
        % are selected among non-migrants in certain years, it is expected
        % to observe unstable coefficient estimates over time.


        % Transfer data to R

    cd('/Users/fgarip/Documents/Clustering Paper 07/Analysis/');

    dlmwrite('coef.txt', coef,',');





% (X) Distribution of Characteristics within and across Clusters %
%     Density Plots - not useful for now %%%%%%%%%%%%%%%%%%%%%%%%%
% 
% cd('/Users/fgarip/Documents/Clustering Paper 07/Analysis/Matlab Analysis/Matlab Plots');
% col = {'b-','r-','g-','y-','m-','k-'};
% 
% for i = 1:size(datac,2) 
%     figure(i);
%     hold on;
%     for j = 1:nocl  
%         [f,xi] = ksdensity(datac(find(clcur==j),i),'kernel','normal');
%         plot(xi,f,col{j});
%     end
%     xlabel(datanamesc{i});
%     ylabel('Density');
%     figure(i);
%     hold off;
%     saveas(gcf,datanamesc{i},'tif'); 
%     %close(i);
% end
%     

% (Y) Visualize clusters using PCA %%%%%%%%%%%%%%%%%
%     Not useful for now 3-dim explain very little %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% [pc, zscores, pcvars] = princomp(datacur, 'econ');
% cumvar = cumsum(pcvars./sum(pcvars)*100)
% pc(:,1:3)
% 
% plot3(zscores(:,1), zscores(:,2), zscores(:,3),'kd')
%     
% col = {'bx','ro','go','ys','mx','kd'};
%     
% figure(1);
% hold on;
%    
% for i = 1:nocl
%     plot3(zscores(find(clcur==i),1), zscores(find(clcur==i),2), zscores(find(clcur==i),3), col{i})
% end
%     
% xlabel('PCA 1');
% ylabel('PCA 2');
% zlabel('PCA 3');
% title('Principal Component Scatter Plot');
% hold off
% 





















