reverse-shooting

Matlab scripts for reverse shooting
Log | Files | Refs | README

Calibrate.m (4240B)

      1 %%% prints values for various values that represent a calibration that's specified in Appendix B1 of the paper. The objective function is based on the transition path. Note that, for each proposed choice of these variables, s_0 and l_0 are found using getells0.m, which is called by solvetransition.m as it gets the transition path.
      2 
      3 %% Initialization
      4 clear all; 
      5 matlabminiscriptspath="/home/[your username]/Documents/ReverseShooting/matlabminiscripts" 
      6 addpath(matlabminiscriptspath)
      7   %%% 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.
      8   %%% Previously: '~/Documents/MATLAB/ChadMatlab/'
      9   %%% addpath just lets matlab know where these files are, much like your $PATH variable in Linux
     10   %%% To be clear, we're still on the same folder, which we can find with the pwd command. 
     11 
     12 %% Files we will work with
     13 if exist('Calibrate.log'); delete('Calibrate.log'); end;
     14 diary Calibrate.log;
     15   %%% The command diary concatenates all subsequent commands and their outputs to whatever is on file Transition.log. If the file doesn't exist, it is created.
     16   %%% This was previously "diary Transition.log;", rather than "Calibrate.log", but was probably a mistake - Nuño.
     17 fprintf(['Calibrate                 ' date]);
     18 disp ' ';
     19 disp ' ';
     20   %%% Some blank space.
     21 
     22 %% Variable definitions
     23 %%% Key values
     24 epsilon=0.4
     25 beta=0.3
     26 gamma=1.5
     27 
     28 phi=5/6
     29 lambda=0.3
     30 
     31 dltabar=1e-4
     32 Nend=1e15
     33 dlta0=5e-4
     34 ubar=1e-2
     35 
     36 %%% Other fixed parameters
     37 rho=.02
     38 alpha=1  %%% 2 percent growth. Nuño: Why does this correspond to 2 percent growth? No idea.
     39 nbar=.01
     40 T=3000
     41 tstep=1
     42 
     43 xguess = [dltabar Nend dlta0 ubar lambda];
     44   %%% This is just an array
     45 
     46 
     47 % Main. Here is where stuff happens
     48 
     49 %% Define options for search function
     50 options = optimoptions('patternsearch', 'UseParallel',true);
     51   %%% options=optimset('Display','on','MaxFunEvals', 100000, 'MaxIter', 10000);
     52 
     53 %% Search a wide space of possibilities. 
     54 xsoln=patternsearch(@(x0) SSRCalibrate(x0, phi, gamma, beta, epsilon, alpha, rho, nbar, T, tstep), xguess, [],[],[],[],[],[],[],options);
     55   %%%xsoln=fminsearch(@(x0) SSRCalibrate(x0, phi, gamma, beta, epsilon, alpha, rho, nbar, T, tstep), xguess, options);
     56   %%% SSRCalibration is defined below (function e=SSRCalibration... is ***
     57   %%% defining the function, counterintuitively
     58   %%% I don't know by what wizardry matlab allows you to define a function
     59   %%% after calling it; this doesn't work with a normal script. 
     60 
     61 %% Consider the possible solutions
     62 xsoln(1)
     63 xsoln(2)
     64 xsoln(3)
     65 xsoln(4)
     66 xsoln(5)
     67 
     68 diary off;
     69   %%% Nuño: Added this last line.
     70 
     71 %% Define our objective function
     72 function e=SSRCalibrate(x0, phi, gamma, beta, epsilon, alpha, rho, nbar, T, tstep)
     73 
     74 dltabar=x0(1); 
     75 
     76 %%% Calculate definitions
     77 %%% x0 is an array, which we take as input, and which will have been
     78   %%% produced by patternsearch.
     79 Nend=x0(2);
     80 dlta0 =x0(3);
     81 ubar =x0(4);
     82 lambda=x0(5);
     83 
     84 [sstar, ellstar, sigmastar, dltastar, ystar, zstar, gs, gc, gh, gdelta] = getsteadystate(dltabar,ubar,epsilon,beta,gamma,dlta0,Nend,alpha,lambda,phi,rho,nbar);
     85   %%% getsteadystate is defined in the function getsteadystate.m
     86 
     87 [t,x,chat,hhat,gdpgrowth,shat,ellhat, dltahat, sigmahat]=solvetransition(dltabar,ubar,epsilon,beta,gamma,dlta0, Nend,alpha,lambda,phi,rho,nbar,T,tstep,0, 0);
     88   %%% solvetransition defined in solvetransition.m
     89   
     90 s=x(:,1);
     91 ell=x(:,2);
     92 sigma=x(:,3);
     93 dlta=x(:,4);
     94 y=x(:,5);
     95 z=x(:,6);
     96 N=x(:,7);
     97 
     98 %%% Calculate distances. 
     99 
    100 ll=ell./(1-ell);
    101 ss=s./(1-s);
    102 
    103 AoverB=(z./y).^(1/(1-phi)).*ss.^(lambda/(1-phi));
    104 coverh=(AoverB.^alpha) .*ll;
    105 c=(((dlta./dltabar).*(coverh.^(-beta))).^(1/(epsilon-beta)))./N;
    106 utilde=ubar.*(c.^(gamma-1))+1/(1-gamma);   % u(c)/u'(c)c
    107 
    108 [minValue,closestIndex] = min(abs(utilde-4));
    109 
    110 m(1)=((utilde(closestIndex) - 4)/4);
    111 m(2)=3*((gdelta-dltahat(1))/gdelta);
    112 m(3)=((ell(closestIndex) -0.95)/0.95);
    113 %m(4)=((N(closestIndex) - 7e9)/7e9)/100;
    114 m(4)=((chat(closestIndex) -0.01)/0.01)/2;
    115 m(5)=5*((sigmahat(closestIndex) - 0.02/0.02));
    116 m(6)=2*((dlta(closestIndex) - 0.001/0.001));
    117 
    118 %%% Return the output of our loss function
    119 m=m';
    120 e=10000*m'*m;  % Sum of squared deviations
    121 
    122 end