%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % 006_markov_chain.m % % % % Max Glonek, University of Adelaide % % 10 June, 2019 % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % This code models graphs as discrete-time Markov chains, and exactly % % calculates absorption probabilities and expected times to absorption % % for all transient states in the chain (unlabelled nodes in the graph), % % for each house and senate in this study. A supplementary list of node % % names/IDs is also generated. Files are saved in the "final data" % % subdirectory. % % % % Taking the 74th House as an example, the following file naming % % conventions are used: % % - absorption probabilities: "H074_prob.csv" % % - expected times to absorption: "H074_time.csv" % % - node names/IDs: "H074_names.csv" % % % % Taking the 110th Senate as an example, the following file naming % % conventions are used: % % - absorption probabilities: "S110_prob.csv" % % - expected times to absorption: "S110_time.csv" % % - node names/IDs: "S110_names.csv" % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % get roles for the house roles = readtable("roles.csv"); roles = roles(find(roles{:,2}=="house"),:); % for each house for i = 74:115 %74 - 115 % read in adjacency matrix and node names (icpsr) adj_raw = readtable(join(["out data/H",num2str(i,"%03.f"),"_mat.csv"],""),"HeaderLines",1); names_raw = readtable(join(["out data/H",num2str(i,"%03.f"),"_names.csv"],"")); adj = adj_raw{:,:}; names = str2double(names_raw.x); % get dems and reps roles_i = roles(find(roles{:,1}==i),:); % get all roles for congress roles_d = roles_i(find(roles_i{:,3}=="D"),:); % get dems for congress icpsr_d = roles_d{:,4}; % get dem icpsrs roles_r = roles_i(find(roles_i{:,3}=="R"),:); % get reps for congress icpsr_r = roles_r{:,4}; % get rep icpsrs % specify new order (unknown, dems, reps) idx_d = find(ismember(names,icpsr_d)); % dem indices idx_r = find(ismember(names,icpsr_r)); % rep indices idx_u = find(~ismember(names,[icpsr_d;icpsr_r])); idx = [idx_u;idx_d;idx_r]; % reorder columns prob = adj(idx,idx); % set lengths u = length(idx_u); % unlabelled d = length(idx_d); % dem r = length(idx_r); % rep l = d + r; % buils transition matrix sums = zeros(u,1); for j = 1:u % get denominators sums(j) = sum(prob(j,:)); end for j = 1:u % row sums 1 for transient states prob(j,:) = prob(j,:)/sums(j); end % fix absorbing states (dem,rep) prob(u+1:u+l,:) = [zeros(l,u),eye(l)]; % check row sums -- good to 13 DP check = zeros(1,u); for j = 1:l check(j) = round(sum(prob(j,:)),13); if check(j) < 1 i end end % Consider form [ R | S ] % [---+---] % [ 0 | I ] % Calculate hitting probabilities and times R = prob(1:u,1:u); S = prob(1:u,u+1:u+l); I = eye(u); P = (I-R)^-1*S; % hitting probabilities T = (I-R)^-1*ones(u,1); % hitting times % calculate class absorption probabilities if l == 2 % if there are only two leaders P2 = P; else if d == 1 % if one dem P_d = P(:,1); % get dem probs P_r = zeros(u,1); % sum rep probs for j = 1:r P_r = P_r + P(:,j+1); end elseif r == 1 % if one rep P_d = zeros(u,1); % sum dem probs for j = 1:d P_d = P_d + P(:,j); end P_r = P(:,d+1); % get rep probs else P_d = zeros(u,1); % sum dem probs for j = 1:d P_d = P_d + P(:,j); end P_r = zeros(u,1); % sum rep probs for j = 1:r P_r = P_r + P(:,j+1); end end P2 = [P_d,P_r]; % join P_d and P_r end % generate unlabelled names names_u = names(idx_u); % save outputs dlmwrite(join(["final data/H",num2str(i,"%03.f"),"_prob.csv"],""),P2,'delimiter',',','precision',16); % write probs dlmwrite(join(["final data/H",num2str(i,"%03.f"),"_time.csv"],""),T,'delimiter',',','precision',16); % write times dlmwrite(join(["final data/H",num2str(i,"%03.f"),"_names.csv"],""),names_u,'delimiter',',','precision',16); % write names end % get roles for the senate roles = readtable("roles.csv"); roles = roles(find(roles{:,2}=="senate"),:); % for each senate for i = 74:115 %74 - 115 % read in adjacency matrix and node names (icpsr) adj_raw = readtable(join(["out data/S",num2str(i,"%03.f"),"_mat.csv"],""),"HeaderLines",1); names_raw = readtable(join(["out data/S",num2str(i,"%03.f"),"_names.csv"],"")); adj = adj_raw{:,:}; names = str2double(names_raw.x); % get dems and reps roles_i = roles(find(roles{:,1}==i),:); % get all roles for congress roles_d = roles_i(find(roles_i{:,3}=="D"),:); % get dems for congress icpsr_d = roles_d{:,4}; % get dem icpsrs roles_r = roles_i(find(roles_i{:,3}=="R"),:); % get reps for congress icpsr_r = roles_r{:,4}; % get rep icpsrs % specify new order (unknown, dems, reps) idx_d = find(ismember(names,icpsr_d)); % dem indices idx_r = find(ismember(names,icpsr_r)); % rep indices idx_u = find(~ismember(names,[icpsr_d;icpsr_r])); idx = [idx_u;idx_d;idx_r]; % reorder columns prob = adj(idx,idx); % set lengths u = length(idx_u); % unlabelled d = length(idx_d); % dem r = length(idx_r); % rep l = d + r; % buils transition matrix sums = zeros(u,1); for j = 1:u % get denominators sums(j) = sum(prob(j,:)); end for j = 1:u % row sums 1 for transient states prob(j,:) = prob(j,:)/sums(j); end % fix absorbing states (dem,rep) prob(u+1:u+l,:) = [zeros(l,u),eye(l)]; % check row sums -- good to 13 DP check = zeros(1,u); for j = 1:l check(j) = round(sum(prob(j,:)),13); if check(j) < 1 i end end % Consider form [ R | S ] % [---+---] % [ 0 | I ] % Calculate hitting probabilities and times R = prob(1:u,1:u); S = prob(1:u,u+1:u+l); I = eye(u); P = (I-R)^-1*S; % hitting probabilities T = (I-R)^-1*ones(u,1); % hitting times % calculate class absorption probabilities if l == 2 % if there are only two leaders P2 = P; else if d == 1 % if one dem P_d = P(:,1); % get dem probs P_r = zeros(u,1); % sum rep probs for j = 1:r P_r = P_r + P(:,j+1); end elseif r == 1 % if one rep P_d = zeros(u,1); % sum dem probs for j = 1:d P_d = P_d + P(:,j); end P_r = P(:,d+1); % get rep probs else P_d = zeros(u,1); % sum dem probs for j = 1:d P_d = P_d + P(:,j); end P_r = zeros(u,1); % sum rep probs for j = 1:r P_r = P_r + P(:,j+1); end end P2 = [P_d,P_r]; % join P_d and P_r end % generate unlabelled names names_u = names(idx_u); % save outputs dlmwrite(join(["final data/S",num2str(i,"%03.f"),"_prob.csv"],""),P2,'delimiter',',','precision',16); % write probs dlmwrite(join(["final data/S",num2str(i,"%03.f"),"_time.csv"],""),T,'delimiter',',','precision',16); % write times dlmwrite(join(["final data/S",num2str(i,"%03.f"),"_names.csv"],""),names_u,'delimiter',',','precision',16); % write names end