TransitoryShock.m (7416B)
1 % TransitoryShock.m 7/28/19 2 % 3 % Show results for a given set of parameter values; show what happens when 4 % growth accelerates (transitory). 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('TransitoryShock.log'); delete('TransitoryShock.log'); end; 21 diary TransitoryShock.log; 22 fprintf(['TransitoryShock ' date]); 23 disp ' '; 24 disp ' '; 25 help TransitoryShock 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, 2); 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, 2); 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 139 idiff = closestIndexShock - ibeforesshock; 140 141 142 % calculate survival probability 143 dltasum = sum(dlta(1:closestIndex))*tstep + dlta(1)/gdelta 144 M = exp(-dltasum) 145 146 147 % calculate survival probability 148 dltasumshock = sum(dlta2(1:closestIndex+idiff))*tstep + dlta2(1)/gdelta 149 Mshock = exp(-dltasumshock) 150 151 152 153 154 155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 156 % Transition plot 157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 158 159 tstart=(yeartoday-600); % Starting point for "nice" plots -- make this time 0 160 t=t-tstart; 161 ii=find(t>=0); 162 163 if (idiff <= 0) 164 ii2 =ii(1:end+idiff); 165 ii=ii(-idiff+1:end); 166 t=t(ii); 167 else 168 ii2=ii(idiff+1:end); 169 ii =ii(1:end-idiff); 170 t=t(ii); 171 end 172 173 174 175 176 177 178 % Note well: t goes in *reverse* order 179 180 % Growth rates 181 figure(1); figsetup; 182 plot(t,y2(ii2),'--','Color',[myblue opacityalt],'LineWidth',lw); hold on; 183 plot(t,y(ii),'--','Color',myblue,'LineWidth',lw); 184 plot(t,z2(ii2),'--','Color',[mygreen opacityalt],'LineWidth',lw); 185 plot(t,z(ii),'--','Color',mygreen,'LineWidth',lw); 186 187 plot(t,chat2(ii2),'-','Color',[myblue opacityalt],'LineWidth',lw); 188 plot(t,chat(ii),'-','Color',myblue, 'LineWidth',lw); 189 190 plot(t,hhat2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw); 191 plot(t,hhat(ii),'-','Color',mygreen,'LineWidth',lw); 192 193 ax=axis; ax(1)=0; ax(2)=t(1); ax(3) = 0;ax(4) = 0.02; axis(ax); 194 chadfig('Time','Growth rate',1,0); 195 makefigwide; 196 set(gca,'YTickLabel',strmat('0 1% 2% 3%')); 197 set(gca,'YTick',[0 .01 .02 .03]); 198 199 text(950,.017,'Safety, $h$','Color',mygreen,'interpreter','latex'); 200 text(1200,.012,'Safety tech, $B$','Color',mygreen,'interpreter','latex'); 201 text(1000,.003,'Consumption, $c$','Color','b','interpreter','latex'); 202 text(1200,.008,'Consumption tech, $A$','Color','b','interpreter','latex'); 203 % text(260,.025,'Per capita GDP','Color',mypurp,'interpreter','latex'); 204 print -depsc ../graphs/TransitoryShockGrowthRates.eps 205 206 207 % Labor allocation 208 figure(2); figsetup; 209 plot(t,100*ell2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw);hold on; 210 plot(t,100*ell(ii),'-','Color',mygreen,'LineWidth',lw); 211 212 plot(t,100*s2(ii2),'-','Color',[mypurp opacityalt],'LineWidth',lw); 213 plot(t,100*s(ii),'-','Color',mypurp,'LineWidth',lw); 214 215 plot(t,100*sigma2(ii2),'-', 'Color',[myblue opacityalt], 'LineWidth',lw); 216 plot(t,100*sigma(ii),'-', 'Color',myblue, 'LineWidth',lw); 217 218 ax=axis; ax(1)=0;ax(2)=t(1);ax(3)=0; ax(4)=100; axis(ax); 219 chadfig('Time','Percent',1,0); 220 makefigwide; 221 222 text(1000,80,mlstring('Share of scientists\\ \ \ in $c$ sector, $s$'),'Color',mypurp,'interpreter','latex'); 223 text(300,75,mlstring('Share of workers\\ \ \ in $c$ sector, $\ell$'),'Color',mygreen,'interpreter','latex'); 224 text(400,15,mlstring('Share of researchers\\ in the population, $\sigma$'),'Color','b','interpreter','latex'); 225 print -depsc ../graphs/TransitoryShockAllocation.eps 226 227 228 % Mortality rate 229 figure(3); figsetup; 230 plot(t,100*dlta2(ii2),'-','Color',[0,0,0, opacityalt],'LineWidth',lw); hold on; 231 plot(t,100*dlta(ii),'-','Color',[0,0,0],'LineWidth',lw); 232 233 ax=axis; ax(1)=0;ax(2)=t(1);ax(3)=0; axis(ax); 234 chadfig('Time','Percent',1,0); 235 makefigwide; 236 % set(gca,'YTickLabel',strmat('0 0.01% 0.02% 0.03% 0.04%')); 237 % set(gca,'YTick',[0 .01 .02 .03 .04]); 238 text(700,.15,'Hazard rate, $\delta$','Color',[0,0,0],'interpreter','latex'); 239 print -depsc ../graphs/TransitoryShockMortality.eps 240 241 242 % Population 243 figure(4); figsetup; 244 plot(t,log(N(ii)),'-','Color',mypurp,'LineWidth',lw); hold on; 245 plot(t,log(N2(ii2)),'-','Color',[mypurp opacityalt],'LineWidth',lw); 246 247 ax=axis; ax(1)=0; axis(ax); 248 chadfig('Time','Log population',1,0); 249 makefigwide; 250 251 252 % Utilde 253 figure(5); figsetup; 254 plot(t,utilde2(ii2),'-','Color',[mygreen opacityalt],'LineWidth',lw); hold on; 255 plot(t,utilde(ii),'-','Color',mygreen,'LineWidth',lw); 256 ax=axis; ax(1)=0; ax(2)=t(1); axis(ax); 257 chadfig('Time','Relative value of life',1,0); 258 makefigwide; 259 print -depsc ../graphs/TransitoryShockValueOfLife.eps 260 261 % % Consumption 262 % figure(5); figsetup; 263 % plot(t,log(c(ii)),'-','Color',myblue,'LineWidth',lw); hold on; 264 % plot(t,log(c2(ii2)),'-','Color',[myblue opacityalt],'LineWidth',lw); 265 % ax=axis; ax(1)=0; axis(ax); 266 % chadfig('Time','Log consumption per capita',1,0); 267 % makefigwide; 268 % 269 % 270 % % Safety 271 % figure(6); figsetup; 272 % plot(t,log(h(ii)),'-','Color',mygreen,'LineWidth',lw); hold on; 273 % plot(t,log(h2(ii2)),'-','Color',[mygreen opacityalt],'LineWidth',lw); 274 % ax=axis; ax(1)=0; axis(ax); 275 % chadfig('Time','Log safety per capita',1,0); 276 % makefigwide; 277 % 278 279 280 281 282 diary off;