Transition.m (6638B)
1 %%% Transition.m 2 %%% Show results for a given set of parameter values. 3 %%% Numerical solution of transition dynamics for the existential risk 4 %%% growth model 5 %%% The system of 6 differential equations 6 %%% x=[s,ell,sigma,dlta,y,z] y==gA z==gB 7 %%% dx=transit1dx(t,x,alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar); 8 9 % Initialization 10 clear all; 11 matlabminiscriptspath="/home/[your username]/Documents/ReverseShooting/matlabminiscripts" 12 addpath(matlabminiscriptspath) 13 %%% The matlabminiscriptspath folder just contains some math utilities which used to be present in Chad Jones' computer, from whose paper Life and Growth (https://web.stanford.edu/~chadj/LifeandGrowthJPE2016.pdf) this code originally comes. 14 %%% Previously: '~/Documents/MATLAB/ChadMatlab/' 15 %%% addpath just lets matlab know where these files are, much like your $PATH variable in Linux 16 %%% To be clear, we're still on the same folder, which we can find with the pwd command. 17 18 %% Files we will work with 19 if exist('Transition.log'); delete('Transition.log'); end; 20 diary Transition.log; 21 fprintf(['Transition ' date]); 22 disp ' '; 23 disp ' '; 24 help Transition 25 26 %% Change font size 27 set(0,'defaultAxesFontSize',13); 28 set(0,'defaultTextFontSize',13); 29 30 %% Graph parameters 31 %%% Colors 32 mygreen=[0 .6 .4]; 33 mypurp=[.8 .1 .6]; 34 myblue=[0 .1 .8]; 35 36 %%% Line width for graphs 37 lw=3; 38 39 %% Variable definitions. 40 %%% Key Values 41 epsilon=0.4 42 beta=0.3 43 gamma=1.5 44 phi=5/6 45 dltabar= 3.8965e-05 46 Nend= 9.2955e+14 47 dlta0= 5.0000e-04 48 ubar= 0.0098 49 lambda= 0.3 50 51 %%% Other fixed parameters 52 rho=.02 53 alpha=1 %%% 2 percent growth. Nuño: Why does this correspond to 2 percent growth? No idea. 54 nbar=.01 55 T=2000 56 tstep=1 57 58 % Main. Here is where stuff happens 59 60 %% Get the steady state 61 [sstar, ellstar, sigmastar, dltastar, ystar, zstar, gs, gc, gh, gdelta] = getsteadystate(dltabar,ubar,epsilon,beta,gamma,dlta0,Nend,alpha,lambda,phi,rho,nbar); 62 63 %% Get the solution 64 ShowResults=1; 65 [t,x,chat,hhat,gdpgrowth,shat,ellhat, dltahat, sigmahat]=solvetransition(dltabar,ubar,epsilon,beta,gamma,dlta0, Nend,alpha,lambda,phi,rho,nbar,T,tstep,ShowResults, 0); 66 67 %% Recover the key variables 68 s=x(:,1); 69 ell=x(:,2); 70 sigma=x(:,3); 71 dlta=x(:,4); 72 y=x(:,5); 73 z=x(:,6); 74 N=x(:,7); 75 76 ll=ell./(1-ell); 77 ss=s./(1-s); 78 79 AoverB=(z./y).^(1/(1-phi)).*ss.^(lambda/(1-phi)); 80 coverh=(AoverB.^alpha) .*ll; 81 c=(((dlta./dltabar).*(coverh.^(-beta))).^(1/(epsilon-beta)))./N; 82 h=c./coverh; 83 utilde=ubar.*(c.^(gamma-1))+1/(1-gamma); % u(c)/u'(c)c 84 uoverv =ll.*beta.*dlta.*utilde+epsilon.*dlta.*utilde; 85 86 gvt = rho - uoverv + dlta; 87 vtilde = utilde./(rho - dlta + gvt); 88 dltavtilde = dlta.*vtilde; 89 90 [minValue,closestIndex] = min(abs(utilde-4)); 91 yeartoday=t(closestIndex) 92 93 %% calculate survival probability 94 dltasum = sum(dlta(1:closestIndex))*tstep + dlta(1)/gdelta 95 Mconditionalontoday = exp(-dltasum) 96 97 98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 99 % Transition plot 100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 101 tstart=yeartoday-600; 102 %%% Starting point for "nice" plots -- make this time 0 103 %%%% Nuño: Another cryptic comment. 104 t=t-tstart; 105 ii=find(t>=0); 106 t=t(ii); 107 ttoday=closestIndex; 108 109 % Note well: t goes in *reverse* order 110 %%% Nuño: What a great idea. 111 112 %% Growth rates 113 figure(1); figsetup; 114 plot(t,y(ii),'--','Color',myblue,'LineWidth',lw); hold on; %blue 115 plot(t,z(ii),'--','Color',mygreen,'LineWidth',lw); 116 plot(t,chat(ii),'-','Color',myblue, 'LineWidth',lw); %blue 117 plot(t,hhat(ii),'-','Color',mygreen,'LineWidth',lw); 118 119 %% Highlight "today" 120 plot(t(ttoday),y(ttoday),'o','Color','b','LineWidth',lw,'MarkerFaceColor','b'); 121 plot(t(ttoday),z(ttoday),'o','Color',mygreen,'LineWidth',lw,'MarkerFaceColor',mygreen); 122 plot(t(ttoday),chat(ttoday),'o','Color','b','LineWidth',lw); 123 plot(t(ttoday),hhat(ttoday),'o','Color',mygreen,'LineWidth',lw); 124 ax=axis; ax(1)=0; ax(2)=t(1); ax(3) = 0; ax(4) = 0.02; axis(ax); 125 chadfig('Time','Growth rate',1,0); 126 makefigwide; 127 set(gca,'YTickLabel',strmat('0 1% 2% 3%')); 128 set(gca,'YTick',[0 .01 .02 .03]); 129 text(950,.017,'Safety, $h$','Color',mygreen,'interpreter','latex'); 130 text(1200,.012,'Safety tech, $B$','Color',mygreen,'interpreter','latex'); 131 text(1000,.003,'Consumption, $c$','Color','b','interpreter','latex'); 132 text(1200,.008,'Consumption tech, $A$','Color','b','interpreter','latex'); 133 %%% text(260,.025,'Per capita GDP','Color',mypurp,'interpreter','latex'); 134 print -depsc ../graphs/TransitionGrowthRates.eps 135 136 137 %% Labor allocation 138 figure(2); figsetup; 139 plot(t,100*ell(ii),'-','Color',mygreen,'LineWidth',lw);hold on; 140 plot(t,100*s(ii),'-','Color',mypurp,'LineWidth',lw); 141 plot(t,100*sigma(ii),'-', 'Color',myblue, 'LineWidth',lw); %blue 142 143 ax=axis; ax(1)=0; ax(2)=t(1);ax(3)=0; ax(4)=100; axis(ax); 144 %%% plot([ttoday ttoday],[0 99],'k--','LineWidth',1); 145 146 %% Highlight today 147 plot(t(ttoday),100*sigma(ttoday),'o','Color',myblue,'LineWidth',lw); 148 plot(t(ttoday),100*ell(ttoday),'o','Color',mygreen,'LineWidth',lw); 149 plot(t(ttoday),100*s(ttoday),'o','Color',mypurp,'LineWidth',lw); 150 chadfig('Time','Percent',1,0); 151 makefigwide; 152 text(1000,80,mlstring('Share of scientists\\ \ \ in $c$ sector, $s$'),'Color',mypurp,'interpreter','latex'); 153 text(300,75,mlstring('Share of workers\\ \ \ in $c$ sector, $\ell$'),'Color',mygreen,'interpreter','latex'); 154 text(400,15,mlstring('Share of researchers\\ in the population, $\sigma$'),'Color','b','interpreter','latex'); 155 print -depsc ../graphs/TransitionAllocation.eps 156 157 158 %% Mortality rate 159 figure(3); figsetup; 160 plot(t,100*dlta(ii),'-','Color',[0,0,0],'LineWidth',lw); hold on; 161 plot(t(ttoday),100*dlta(ttoday),'o','Color',[0,0,0],'LineWidth',lw); 162 ax=axis; ax(1)=0;ax(2)=t(1);ax(3)=0; axis(ax); 163 chadfig('Time','Percent',1,0); 164 makefigwide; 165 %%% set(gca,'YTickLabel',strmat('0 0.01% 0.02% 0.03% 0.04%')); 166 %%% set(gca,'YTick',[0 .01 .02 .03 .04]); 167 text(700,.15,'Hazard rate, $\delta$','Color',[0,0,0],'interpreter','latex'); 168 print -depsc ../graphs/TransitionMortality.eps 169 170 171 %% Population 172 figure(4); figsetup; 173 plot(t,log(N(ii)),'-','Color',mypurp,'LineWidth',lw); hold on; 174 %%% plot(t,log(N2(ii2)),'-','Color',[mypurp 0.5],'LineWidth',lw); 175 176 ax=axis; ax(1)=0; axis(ax); 177 chadfig('Time','Log population',1,0); 178 makefigwide; 179 180 %% Consumption 181 figure(5); figsetup; 182 plot(t,log(c(ii)),'-','Color',myblue,'LineWidth',lw); hold on; 183 %%% plot(t,log(c2(ii2)),'-','Color',[myblue 0.5],'LineWidth',lw); 184 ax=axis; ax(1)=0; axis(ax); 185 chadfig('Time','Log consumption per capita',1,0); 186 makefigwide; 187 188 %% Safety 189 figure(6); figsetup; 190 plot(t,log(h(ii)),'-','Color',mygreen,'LineWidth',lw); hold on; 191 %%% plot(t,log(h2(ii2)),'-','Color',[mygreen 0.5],'LineWidth',lw); 192 ax=axis; ax(1)=0; axis(ax); 193 chadfig('Time','Log safety per capita',1,0); 194 makefigwide; 195 196 diary off;