Acceleration.m (7202B)
1 % Acceleration.m 7/28/19 2 % 3 % Show results for a given set of parameter values; show what happens when 4 % growth accelerates. 5 % 6 % Numerical solution of transition dynamics for the existential risk 7 % growth model 8 % 9 % The system of 6 differential equations 10 % 11 % x=[s,ell,sigma,dlta,y,z] y==gA z==gB 12 % dx=transit1dx(t,x,alpha,epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar); 13 % 14 15 clear all; %startup 16 17 addpath('~/Documents/MATLAB/ChadMatlab/') 18 19 20 if exist('Acceleration.log'); delete('Acceleration.log'); end; 21 diary Acceleration.log; 22 fprintf(['Acceleration ' date]); 23 disp ' '; 24 disp ' '; 25 help Acceleration 26 27 % Change font size 28 set(0,'defaultAxesFontSize',13); 29 set(0,'defaultTextFontSize',13); 30 31 % Parameters 32 mygreen=[0 .6 .4]; 33 mypurp=[.8 .1 .6]; 34 myblue=[0 .1 .8]; 35 36 lw=2; 37 opacityalt=0.4; 38 39 40 % Key Values 41 epsilon=0.4 42 beta=0.3 43 gamma=1.5 44 45 phi=5/6 46 47 dltabar= 3.8965e-05 48 Nend= 9.2955e+14 49 dlta0= 5.0000e-04 50 ubar= 0.0098 51 lambda= 0.3 52 53 54 55 % Other fixed parameters 56 rho=.02 57 alpha=1 % 2 percent growth 58 nbar=.01 59 T=2000 60 tstep=1 61 62 [sstar, ellstar, sigmastar, dltastar, ystar, zstar, gs, gc, gh, gdelta] = getsteadystate(dltabar,ubar,epsilon,beta,gamma,dlta0,Nend,alpha,lambda,phi,rho,nbar); 63 64 % Get solution 65 ShowResults=1; 66 [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); 67 68 69 % Recover the key variables 70 s=x(:,1); 71 ell=x(:,2); 72 sigma=x(:,3); 73 dlta=x(:,4); 74 y=x(:,5); 75 z=x(:,6); 76 N=x(:,7); 77 78 79 ll=ell./(1-ell); 80 ss=s./(1-s); 81 82 AoverB=(z./y).^(1/(1-phi)).*ss.^(lambda/(1-phi)); 83 coverh=(AoverB.^alpha) .*ll; 84 c=(((dlta./dltabar).*(coverh.^(-beta))).^(1/(epsilon-beta)))./N; 85 h=c./coverh; 86 utilde=ubar.*(c.^(gamma-1))+1/(1-gamma); % u(c)/u'(c)c 87 88 89 [minValue,closestIndex] = min(abs(utilde-4)) 90 yeartoday=t(closestIndex) 91 92 93 ibeforesshock = 1300; 94 95 %parameters to calibrate 96 scal = s(ibeforesshock); 97 ellcal = ell(ibeforesshock); 98 sigmacal = sigma(ibeforesshock); 99 dltacal = dlta(ibeforesshock); 100 ycal = y(ibeforesshock); 101 zcal = z(ibeforesshock); 102 Ncal = N(ibeforesshock); 103 104 105 xguess = [dlta0 Nend]; 106 107 calibratedparam=FindAlternativePath(xguess, alpha, epsilon,beta,lambda,phi,dltabar,ubar,gamma,rho,nbar, Ncal, scal, ellcal, sigmacal, dltacal, ycal, zcal, T, tstep, 1); 108 109 dltanew = calibratedparam(1) 110 Nnew= calibratedparam(2) 111 112 113 % Get solution for shock 114 [t2,x2,chat2,hhat2,gdpgrowth2,shat2,ellhat2, dltahat2, sigmahat2]=solvetransition(dltabar,ubar,epsilon,beta,gamma,dltanew,Nnew,alpha,lambda,phi,rho,nbar,T,tstep,ShowResults, 1); 115 116 % Recover the key variables 117 s2=x2(:,1); 118 ell2=x2(:,2); 119 sigma2=x2(:,3); 120 dlta2=x2(:,4); 121 y2=x2(:,5); 122 z2=x2(:,6); 123 N2=x2(:,7); 124 125 [minValue,closestIndexShock] = min(abs(dlta2-dltacal)); 126 127 128 ll2=ell2./(1-ell2); 129 ss2=s2./(1-s2); 130 131 132 AoverB2=(z2./y2).^(1/(1-phi)).*ss2.^(lambda/(1-phi)); 133 coverh2=(AoverB2.^alpha) .*ll2; 134 c2=(((dlta2./dltabar).*(coverh2.^(-beta))).^(1/(epsilon-beta)))./N2; 135 h2=c2./coverh2; 136 utilde2=ubar.*(c2.^(gamma-1))+1/(1-gamma); % u(c)/u'(c)c 137 138 idiff = closestIndexShock - ibeforesshock; 139 140 141 % calculate survival probability 142 dltasum = sum(dlta(1:closestIndex))*tstep + dlta(1)/gdelta 143 M = exp(-dltasum) 144 145 146 % calculate survival probability 147 dltasumshock = sum(dlta2(1:closestIndex+idiff))*tstep + dlta2(1)/gdelta 148 Mshock = exp(-dltasumshock) 149 150 151 152 153 154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 155 % Transition plot 156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 157 158 tstart=(yeartoday-600); % Starting point for "nice" plots -- make this time 0 159 t=t-tstart; 160 ii=find(t>=0); 161 162 ii2 =ii(1:end+idiff); 163 ii=ii(-idiff+1:end); 164 t=t(ii); 165 166 167 168 169 % Note well: t goes in *reverse* order 170 171 % Growth rates 172 figure(1); figsetup; 173 plot(t,y2(ii2),'--','Color',[myblue opacityalt],'LineWidth',lw); hold on; 174 plot(t,y(ii),'--','Color',myblue,'LineWidth',lw); 175 plot(t,z2(ii2),'--','Color',[mygreen opacityalt],'LineWidth',lw); 176 plot(t,z(ii),'--','Color',mygreen,'LineWidth',lw); 177 178 plot(t,chat2(ii2),'-','Color',[myblue opacityalt],'LineWidth',lw); 179 plot(t,chat(ii),'-','Color',myblue, 'LineWidth',lw); 180 181 plot(t,hhat2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw); 182 plot(t,hhat(ii),'-','Color',mygreen,'LineWidth',lw); 183 184 ax=axis; ax(1)=0; ax(2)=t(1); ax(3) = 0;ax(4) = 0.02; axis(ax); 185 chadfig('Time','Growth rate',1,0); 186 makefigwide; 187 set(gca,'YTickLabel',strmat('0 1% 2% 3%')); 188 set(gca,'YTick',[0 .01 .02 .03]); 189 text(950,.017,'Safety, $h$','Color',mygreen,'interpreter','latex'); 190 text(1200,.012,'Safety tech, $B$','Color',mygreen,'interpreter','latex'); 191 text(1000,.003,'Consumption, $c$','Color','b','interpreter','latex'); 192 text(1200,.008,'Consumption tech, $A$','Color','b','interpreter','latex'); 193 % text(260,.025,'Per capita GDP','Color',mypurp,'interpreter','latex'); 194 print -depsc ../graphs/AccelerationGrowthRates.eps 195 196 197 % Labor allocation 198 figure(2); figsetup; 199 plot(t,100*ell2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw);hold on; 200 plot(t,100*ell(ii),'-','Color',mygreen,'LineWidth',lw); 201 202 plot(t,100*s2(ii2),'-','Color',[mypurp opacityalt],'LineWidth',lw); 203 plot(t,100*s(ii),'-','Color',mypurp,'LineWidth',lw); 204 205 plot(t,100*sigma2(ii2),'-', 'Color',[myblue opacityalt], 'LineWidth',lw); %blue 206 plot(t,100*sigma(ii),'-', 'Color',myblue, 'LineWidth',lw); %blue 207 208 ax=axis; ax(1)=0;ax(2)=t(1);ax(3)=0; ax(4)=100; axis(ax); 209 chadfig('Time','Percent',1,0); 210 makefigwide; 211 text(1000,80,mlstring('Share of scientists\\ \ \ in $c$ sector, $s$'),'Color',mypurp,'interpreter','latex'); 212 text(300,75,mlstring('Share of workers\\ \ \ in $c$ sector, $\ell$'),'Color',mygreen,'interpreter','latex'); 213 text(400,15,mlstring('Share of researchers\\ in the population, $\sigma$'),'Color','b','interpreter','latex'); 214 print -depsc ../graphs/AccelerationAllocation.eps 215 216 217 % Mortality rate 218 figure(3); figsetup; 219 plot(t,100*dlta2(ii2),'-','Color',[0,0,0, opacityalt],'LineWidth',lw); hold on; 220 plot(t,100*dlta(ii),'-','Color',[0,0,0],'LineWidth',lw); 221 222 ax=axis; ax(1)=0;ax(2)=t(1);ax(3)=0; axis(ax); 223 chadfig('Time','Percent',1,0); 224 makefigwide; 225 % set(gca,'YTickLabel',strmat('0 0.01% 0.02% 0.03% 0.04%')); 226 % set(gca,'YTick',[0 .01 .02 .03 .04]); 227 text(700,.15,'Hazard rate, $\delta$','Color',[0,0,0],'interpreter','latex'); 228 print -depsc ../graphs/AccelerationMortality.eps 229 230 231 232 % % Mortality rate 233 % figure(3); figsetup; 234 % plot(t,100*dlta2(ii2),'-','Color',[0,0,0, opacityalt],'LineWidth',lw); hold on; 235 % plot(t,100*dlta(ii),'-','Color',[0,0,0],'LineWidth',lw); 236 % 237 % ax=axis; ax(1)=550; ax(2) = 700; ax(3)=0; axis(ax); 238 % chadfig('Time','Mortality rate in percent',1,0); 239 % makefigwide; 240 % % set(gca,'YTickLabel',strmat('0 0.01% 0.02% 0.03% 0.04%')); 241 % % set(gca,'YTick',[0 .01 .02 .03 .04]); 242 % print -depsc ../graphs/AccelerationMortalityZoomedIn.eps 243 244 245 % Population 246 figure(4); figsetup; 247 plot(t,log(N(ii)),'-','Color',mypurp,'LineWidth',lw); hold on; 248 plot(t,log(N2(ii2)),'-','Color',[mypurp opacityalt],'LineWidth',lw); 249 250 ax=axis; ax(1)=0; axis(ax); 251 chadfig('Time','Log population',1,0); 252 makefigwide; 253 254 % Utilde 255 figure(5); figsetup; 256 plot(t,utilde2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw);hold on; 257 plot(t,utilde(ii),'-','Color',mygreen,'LineWidth',lw); 258 ax=axis; ax(1)=0;ax(2)=t(1); axis(ax); 259 chadfig('Time','Relative value of life',1,0); 260 makefigwide; 261 print -depsc ../graphs/AccelerationValueOfLife.eps 262 263 264 265 266 267 diary off;