% L Q O C _ U S . M % ******************* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Linear Quadratic Optimal Control of A Small Model of the U.S. Economy : % % ********************************************************************** % % % % x*_t - x*_t-1 = A x*_t-1 + B u*_t + C z_t % % x*_0 = xo % % J* --> min % % % % Author: Osvald Vasicek % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all load us_state % Historical State Data x= XDAT; % Historical Control Data u= UDAT; % Historical Exogenous Variable Data zz= ZDAT; % Data Range Ndat=length(x); n=19; % Number of States r=4; % Number of Controls s=4; % Number of Exogenous Variables % Upper and Lower Limit for Optimalization Nf=25; %<---<<< Ndat, Ndat-10 % Choice of Upper Limit for Optimalization ! ! ! Ns=5; %<---<<< 100; %<---<<< 5, 100 Ndat-48, Ndat -151 % Choice of Lower Limit for Optimalization ! ! ! disp(' '); disp(' '); disp(' Start and Finish of Optimal Control'); disp(' ***********************************'); disp([' t_s = ',num2str(t(Ns))]); disp([' t_f = ',num2str(t(Nf))]); % Historical Value x_hist=x(Ns:Nf,:)'; u_hist=u(Ns:Nf,:)'; N=max(size(x_hist)); % Regress Aproximation of Variable Development through Exponential Functions % -------------------------------------------------------------------------- t_r=t(Ns-4:Nf)'; t4=(1:length(t_r))'; Ar=[ones(size((Ns-4:Nf)')),t4]; [par1,i_par1,e1,i_e1,Stat1]=regress(log(YD(Ns-4:Nf)),Ar,0.05); YDr(Ns-4:Nf)=exp(log(YD(Ns-4:Nf))-e1); [par2,i_par2,e2,i_e2,Stat2]=regress(log(CC(Ns-4:Nf)),Ar,0.05); Cr(Ns-4:Nf)=exp(log(CC(Ns-4:Nf))-e2); [par3,i_par3,e3,i_e3,Stat3]=regress(log(TAX(Ns-4:Nf)),Ar,0.05); TAXr(Ns-4:Nf)=exp(log(TAX(Ns-4:Nf))-e3); [par4,i_par4,e4,i_e4,Stat4]=regress(log(CC(Ns-4:Nf)),Ar,0.05); Cr(Ns-4:Nf)=exp(log(CC(Ns-4:Nf))-e4); [par5,i_par5,e5,i_e5,Stat5]=regress(log(INR(Ns-4:Nf)),Ar,0.05); INRr(Ns-4:Nf)=exp(log(INR(Ns-4:Nf))-e5); [par6,i_par6,e6,i_e6,Stat6]=regress(log(IR(Ns-4:Nf)),Ar,0.05); IRr(Ns-4:Nf)=exp(log(IR(Ns-4:Nf))-e6); [par7,i_par7,e7,i_e7,Stat7]=regress(IIN(Ns-4:Nf),Ar,0.05); IINr(Ns-4:Nf)=IIN(Ns-4:Nf)-e7; [par8,i_par8,e8,i_e8,Stat8]=regress(log(RS(Ns-4:Nf)),Ar,0.05); RSr(Ns-4:Nf)=exp(log(RS(Ns-4:Nf))-e8); [par9,i_par9,e9,i_e9,Stat9]=regress(log(RL(Ns-4:Nf)),Ar,0.05); RLr(Ns-4:Nf)=exp(log(RL(Ns-4:Nf))-e9); [par10,i_par10,e10,i_e10,Stat10]=regress(INF(Ns-4:Nf),Ar,0.05); INFr(Ns-4:Nf)=INF(Ns-4:Nf)-e10; [par11,i_par11,e11,i_e11,Stat11]=regress(WINF(Ns-4:Nf),Ar,0.05); WINFr(Ns-4:Nf)=WINF(Ns-4:Nf)-e11; [par12,i_par12,e12,i_e12,Stat12]=regress(log(UR(Ns-4:Nf)),Ar,0.05); URr(Ns-4:Nf)=exp(log(UR(Ns-4:Nf))-e12); [par13,i_par13,e13,i_e13,Stat13]=regress(log(G(Ns-4:Nf)),Ar,0.05); Gr(Ns-4:Nf)=exp(log(G(Ns-4:Nf))-e13); [par14,i_par14,e14,i_e14,Stat14]=regress(RGM(Ns-4:Nf),Ar,0.05); RGMr(Ns-4:Nf)=RGM(Ns-4:Nf)-e14; [par15,i_par15,e15,i_e15,Stat15]=regress(log(TR(Ns-4:Nf)),Ar,0.05); TRr(Ns-4:Nf)=exp(log(TR(Ns-4:Nf))-e15); [par16,i_par16,e16,i_e16,Stat16]=regress(log(WLTH(Ns-4:Nf)),Ar,0.05); WLTHr(Ns-4:Nf)=exp(log(WLTH(Ns-4:Nf))-e16); % S e t t i n g o f N o m i n a l T r a j e c t o r i e s % ----------------------------------------------------------------------- % Nominal Trajectories: Actual Development: Aproximation Development: % _______________________________________________________________________________ x_nom(1,:) = GNPP(Ns:Nf)'; %...x(Ns:Nf,1)'; %...GNPP(Ns:Nf)'; % 1. GNPP x_nom(2,:) = YDr(Ns:Nf); %...x(Ns:Nf,2)'; %...YDr(Ns:Nf); % 2. YD x_nom(3,:) = TAXr(Ns:Nf); %...x(Ns:Nf,3)'; %...TAXr(Ns:Nf); % 3. TAX x_nom(4,:) = Cr(Ns:Nf); %...x(Ns:Nf,4)'; %...Cr(Ns:Nf); % 4. C x_nom(5,:) = INRr(Ns:Nf); %...x(Ns:Nf,5)'; %...INRr(Ns:Nf); % 5. INR x_nom(6,:) = IRr(Ns:Nf); %...x(Ns:Nf,6)'; %...IRr(Ns:Nf); % 6. IR x_nom(7,:) = IINr(Ns:Nf); %...x(Ns:Nf,7)'; %...IINr(Ns:Nf); % 7. IIN x_nom(8,:) = RSr(Ns:Nf); %...x(Ns:Nf,8)'; %...RSr(Ns:Nf); % 8. RS x_nom(9,:) = RLr(Ns:Nf); %...x(Ns:Nf,9)'; %...RLr(Ns:Nf); % 9. RL x_nom(10,:)= INFr(Ns:Nf); %...x(Ns:Nf,10)'; %...INFr(Ns:Nf); % 10.INF x_nom(11,:)= WINFr(Ns:Nf); %...x(Ns:Nf,11)'; %...WINFr(Ns:Nf); % 11.WINF x_nom(12,:)= URr(Ns:Nf); %...x(Ns:Nf,12)'; %...URr(Ns:Nf); % 12.UR x_nom(13,:)= RLr(Ns-1:Nf-1); %...x(Ns:Nf,13)'; %...RLr(Ns-1:Nf-1); % 13.RL_1 x_nom(14,:)= INFr(Ns-1:Nf-1); %...x(Ns:Nf,14)'; %...INFr(Ns-1:Nf-1); % 14.INF_1 x_nom(15,:)= URr(Ns-1:Nf-1); %...x(Ns:Nf,15)'; %...URr(Ns-1:Nf-1); % 15.UR_1 x_nom(16,:)= RLr(Ns-2:Nf-2); %...x(Ns:Nf,16)'; %...RLr(Ns-2:Nf-2); % 16.RL_2 x_nom(17,:)= INFr(Ns-2:Nf-2); %...x(Ns:Nf,17)'; %...INFr(Ns-2:Nf-2); % 16.INF_2 x_nom(18,:)= URr(Ns-2:Nf-2); %...x(Ns:Nf,18)'; %...URr(Ns-2:Nf-2); % 17.UR_1 x_nom(19,:)= RLr(Ns-3:Nf-3); %...x(Ns:Nf,19)'; %...RLr(Ns-3:Nf-3); % 19.RL_3 u_nom(1,:)= u(Ns:Nf,1)'; %...u(Ns:Nf,1)'; %...Gr(Ns:Nf); u_nom(2,:)= u(Ns:Nf,2)'; %...u(Ns:Nf,2)'; %...RGMr(Ns:Nf); u_nom(3,:)= u(Ns:Nf,3)'; %...u(Ns:Nf,3)'; %...TRr(Ns:Nf); u_nom(4,:)= u(Ns:Nf,4)'; %...u(Ns:Nf,4)'; %...WLTHr(Ns:Nf); z=zz(Ns:Nf,:)'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ATTENTION !!! % % ############# % % If we choose historical/aproximational Development for a "nominal" Devel-% % opment of the given State, we also have to choose the same lagged Devel- % % opment for A l l S t a t e s created by the lag of this State. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Simulation of State Development without Control % ---------------------------------------------- % Initial Condition x_sim_0=x_nom(:,1); x_sim(1:n,1)=x_sim_0; for j = 1:N-1, disp(' d_x_sim_i, x_sim_i'); disp(' i ='); disp(j) d_x_sim(:,j+1) = A*x_sim(:,j) + B*u_hist(:,j+1) + C*z(:,j+1); x_sim(:,j+1) = x_sim(:,j) + d_x_sim(:,j+1); end Q=zeros(19,19); % M a k r o e c o n o m i c S t r a t e g y f o r C o n t r o l % ----------------------------------------------------------------------- % Election of Weight for States % 1 2 3 4 5 6 7 8 9 10 11 12 % GNP YD TAX C INR IR IIN RS RL INF WINF UR % _______________________________________________________________________ w_Q_0=[ 1e-6 0 0 0 0 0 0 0 0 0 0 0 ]; w_Q_1=[ 5e+6 0 0 0 0 0 0 1e+2 0 5e+10 0 0 ]; w_Q_2=[ 1e+5 0 0 0 0 0 0 1e+2 0 1e+3 0 1e+8 ]; % Pindyck "run 1" w_Q_3=[ 0 0 0 1 6 15 0 0 0 6 0 4e+6 ]; % Pindyck "run 2" w_Q_4=[ 100 0 0 1 6 150 0 0 0 120 0 4e+6 ]; w_Q_5=[ 1e+3 0 0 1 6 150 0 1e+3 0 120 0 4e+6 ]; % Choice of Weights for Control % G RGM TR WLTH % _______________________ w_R_0=[ 1e+5 1e+5 1e+5 1e+5]; w_R_1=[ 1e+3 1e+4 8e+6 8e+5]; w_R_2=[ 1e+3 1e+4 1e+5 5e+5]; % Pindyck "run 1" w_R_3=[ 3 300 30 5e+5]; % Pindyck "run 1" w_R_4=[ 3 15 60 5e+5]; w_R_5=[ 3 60 30 5e+5]; % Strategy decision w_R = w_R_2; w_Q = w_Q_2; format bank disp(' '); disp(' '); disp(' Diagonal of Matrix Q'); disp(' GNP YD TAX C INR IR '); disp([num2str(w_Q(1,1:6))]) disp(' IIN RS RL INF WINF UR '); disp([num2str(w_Q(1,7:12))]) disp(' '); disp(' Diagonal of Matrix R'); disp('G RGM TR WLTH '); disp([num2str(w_R)]) format bank % Matrix Q Q(1:12,1:12)=diag(w_Q); % Matrix R R=diag(w_R); format % O P T I M A L C O N T R O L % ############################################## % Boundary Conditions x_opt_0=x_nom(:,1); x_opt_N=x_nom(:,N); K_N = Q; g_N = -Q*x_nom(:,N); % (i) Solving of Riccati Equation for K_i % ************************************* Io =eye(n); K=zeros(n,N*n); % End Conditions K(1:n,(N-1)*n+1:N*n)= K_N; g(1:n,N) = g_N; % Solution K_i for j=N-1:-1:1, disp(' K_i, g_i'); disp(' i ='); disp(j) K1 = K(1:n,(j*n+1):(j+1)*n); K(1:n,((j-1)*n+1):j*n)=Q+(Io+A)'*(K1-K1*B*inv(R+B'*K1*B)*B'*K1)*(Io+A); % (ii) Solving of Determination Equation for g_i % ***************************************** g1 = g(1:n,j+1); g(1:n,j)= -(Io+A)'*(K1-K1*B*inv(R+B'*K1*B)*B'*K1)*B*inv(R)*B'*g1... +(Io+A)'*g1... +(Io+A)'*(K1-K1*B*inv(R+B'*K1*B)*B'*K1)*(B*u_nom(1:r,j+1)+C*z(1:s,j+1))... -Q*x_nom(1:n,j); end % (iii) Solving of Optimal Control u_opt_i and Optimal State x_opt_i % ************************************************************** % Initial Condition x_opt(1:n,1)=x_opt_0; for j = 1:N-1, disp(' u_opt_i, x_opt_i'); disp(' i ='); disp(j) K1 = K(1:n,((j-1)*n+1):(j)*n); g1 = g(1:n,j+1); % Optimal Control u_opt_i u_opt(1:r,j+1)=-inv(R+B'*K1*B)*B'*K1*(Io+A)*x_opt(1:n,j)... +inv(R+B'*K1*B)*B'*K1*B*inv(R)*B'*g1... -inv(R)*B'*g1-inv(R+B'*K1*B)*B'*K1... *(B*u_nom(1:r,j+1)+C*z(1:s,j+1))+u_nom(1:r,j+1); % Optimal Difference of State d_x_opt_i+1 = x_opt_i+1 and State x_opt_i d_x_opt(1:n,j+1) = A*x_opt(1:n,j) + B*u_opt(1:r,j+1) + C*z(1:s,j+1); x_opt(1:n,j+1) = x_opt(1:n,j)+ d_x_opt(1:n,j+1); end u_op = u_opt'; u_nm = u_nom'; x_op = x_opt'; x_sm = x_sim'; x_nm = x_nom'; set(0,'DefaultAxesFontSize',8); set(0,'DefaultLineLineWidth',1); figure(1) clf set(gcf,'Name',' Model States x1(t), ..., x12(t)'); to=t(Ns:Nf); [Nu,n2]=size(u_op); subplot(4,3,1) plot(to,x_nm(:,1),'b.-'); hold on plot(to(2:N),x_op(2:N,1),'c*-'); hold on plot(to,GNP(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,1),'g-.'); hold off legend('nom','opt','orig','sim') axis([min(to)-0.25 max(to)+0.25... min([min(GNP(Ns:Nf)),min(x_sm(:,1)),min(x_op(:,1)),min(x_nm(:,1))])-100.,... max([max(GNP(Ns:Nf)),max(x_sm(:,1)),max(x_op(:,1)),max(x_nm(:,1))])+100]); title('Gross National Product GNP'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,2) plot(to,x_nm(:,2),'b.-'); hold on plot(to(2:N),x_op(2:N,2),'c*-'); hold on plot(to,YD(Ns:Nf),'r-'); hold on plot(to,x_sm(:,2),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(YD(Ns:Nf)),min(x_sm(:,2)),min(x_op(:,2)),min(x_nm(:,2))])-100.,... max([max(YD(Ns:Nf)),max(x_sm(:,2)),max(x_op(:,2)),max(x_nm(:,2))])+100]); title('Disposable Income YD'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,3) plot(to,x_nm(:,3),'b.-'); hold on plot(to(2:N),x_op(2:N,3),'c*-'); hold on plot(to,TAX(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,3),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(TAX(Ns:Nf)),min(x_sm(:,3)),min(x_op(:,3)),min(x_nm(:,3))])-50.,... max([max(TAX(Ns:Nf)),max(x_sm(:,3)),max(x_op(:,3)),max(x_nm(:,3))])+50]); title('Tax Revenue TAX'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,4) plot(to,x_nm(:,4),'b.-'); hold on plot(to(2:N),x_op(2:N,4),'c*-'); hold on plot(to,CC(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,4),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(CC(Ns:Nf)),min(x_sm(:,4)),min(x_op(:,4)),min(x_nm(:,4))])-100.,... max([max(CC(Ns:Nf)),max(x_sm(:,4)),max(x_op(:,4)),max(x_nm(:,4))])+100]); title('Consumption C'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,5) plot(to,x_nm(:,5),'b.-'); hold on plot(to(2:N),x_op(2:N,5),'c*-'); hold on plot(to,INR(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,5),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(INR(Ns:Nf)),min(x_sm(:,5)),min(x_op(:,5)),min(x_nm(:,5))])-25.,... max([max(INR(Ns:Nf)),max(x_sm(:,5)),max(x_op(:,5)),max(x_nm(:,5))])+25]); title('Nonresidential Investment INR'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,6) plot(to,x_nm(:,6),'b.-'); hold on plot(to(2:N),x_op(2:N,6),'c*-'); hold on plot(to,IR(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,6),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(IR(Ns:Nf)),min(x_sm(:,6)),min(x_op(:,6)),min(x_nm(:,6))])-15.,... max([max(IR(Ns:Nf)),max(x_sm(:,6)),max(x_op(:,6)),max(x_nm(:,6))])+15]); title('Residential Investment IR'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,7) plot(to,x_nm(:,7),'b.-'); hold on plot(to(2:N),x_op(2:N,7),'c*-'); hold on plot(to,IIN(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,7),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(IIN(Ns:Nf)),min(x_sm(:,7)),min(x_op(:,7)),min(x_nm(:,7))])-25.,... max([max(IIN(Ns:Nf)),max(x_sm(:,7)),max(x_op(:,7)),max(x_nm(:,7))])+25]); title('Inventory Investment IIN'); ylabel('[ billions of 1982 USD ]'); grid subplot(4,3,8) plot(to,x_nm(:,8),'b.-'); hold on plot(to(2:N),x_op(2:N,8),'c*-'); hold on plot(to,RS(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,8),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(RS(Ns:Nf)),min(x_sm(:,8)),min(x_op(:,8)),min(x_nm(:,8))])-2.,... max([max(RS(Ns:Nf)),max(x_sm(:,8)),max(x_op(:,8)),max(x_nm(:,8))])+2.]); title('Short-term Interest Rate RS'); ylabel('[ pecent per year ]'); grid subplot(4,3,9) plot(to,x_nm(:,9),'b.-'); hold on plot(to(2:N),x_op(2:N,9),'c*-'); hold on plot(to,RL(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,9),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(RL(Ns:Nf)),min(x_sm(:,9)),min(x_op(:,9)),min(x_nm(:,9))])-1.25,... max([max(RL(Ns:Nf)),max(x_sm(:,9)),max(x_op(:,9)),max(x_nm(:,9))])+1.25]); title('Long-term Interest Rate RL'); ylabel('[ pecent per year]'); grid subplot(4,3,10) plot(to,x_nm(:,10),'b.-'); hold on plot(to(2:N),x_op(2:N,10),'c*-'); hold on plot(to,INF(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,10),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(INF(Ns:Nf)),min(x_sm(:,10)),min(x_op(:,10)),min(x_nm(:,10))])-1.,... max([max(INF(Ns:Nf)),max(x_sm(:,10)),max(x_op(:,10)),max(x_nm(:,10))])+1]); title('Inflation Rate INF'); ylabel('[ pecent per year ]'); grid subplot(4,3,11) plot(to,x_nm(:,11),'b.-'); hold on plot(to(2:N),x_op(2:N,11),'c*-'); hold on plot(to,WINF(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,11),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(WINF(Ns:Nf)),min(x_sm(:,11)),min(x_op(:,11)),min(x_nm(:,11))])-1.,... max([max(WINF(Ns:Nf)),max(x_sm(:,11)),max(x_op(:,11)),max(x_nm(:,11))])+1]); title('Rate of Growth of Wage Rate WINF'); ylabel('[ pecent per year ]'); grid subplot(4,3,12) plot(to,x_nm(:,12),'b.-'); hold on plot(to,x_op(:,12),'c*-'); hold on plot(to,UR(Ns:Nf,:)','r-'); hold on plot(to,x_sm(:,12),'g-.'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(UR(Ns:Nf)),min(x_sm(:,12)),min(x_op(:,12)),min(x_nm(:,12))])-0.5,... max([max(UR(Ns:Nf)),max(x_sm(:,12)),max(x_op(:,12)),max(x_nm(:,12))])+0.5]); title('Unemployment Rate UR'); ylabel('[ pecent ]'); grid [Nu,n2]=size(u_op); figure(2) clf set(gcf,'Name','Control and Exogenous Variables i.e. u1(t), ..., u4(t) end z1(t), ..., z3(t) '); subplot(3,3,1) plot(to,u_nm(:,1),'b.-'); hold on plot(to(2:Nu),u_op(2:Nu,1),'c*-'); hold on plot(to,G(Ns:Nf),'r-'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(G(Ns:Nf)),min(u_nm(:,1)),min(u_op(2:Nu,1))])-7,... max([max(G(Ns:Nf)),max(u_nm(:,1)),max(u_op(2:Nu,1))])+7]); title('Government Spending G'); ylabel('[ billions of 1982 USD ]'); grid subplot(3,3,2) plot(to,u_nm(:,2),'b.-'); hold on plot(to(2:Nu),u_op(2:Nu,2),'c*-'); hold on plot(to,RGM(Ns:Nf),'r-'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(RGM(Ns:Nf)),min(u_nm(:,2)),min(u_op(2:Nu,2))])-4,... max([max(RGM(Ns:Nf)),max(u_nm(:,2)),max(u_op(2:Nu,2))])+4]); title('Growth Rate of Money Supply RGM'); ylabel('[ pecent per year ]'); grid subplot(3,3,3) plot(to,u_nm(:,3),'b.-'); hold on plot(to(2:Nu),u_op(2:Nu,3),'c*-'); hold on plot(to,TR(Ns:Nf),'r-'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(TR(Nf:Ns)),min(u_nm(:,3)),min(u_op(2:Nu,3))])-5.0,... max([max(TR(Nf:Ns)),max(u_nm(:,3)),max(u_op(2:Nu,3))])+5.0]); title('Transfer Payments TR'); ylabel('[ billions of 1982 USD ]'); grid subplot(3,3,4) plot(to,u_nm(:,4),'b.-'); hold on plot(to(2:Nu),u_op(2:Nu,4),'c*-'); hold on plot(to,WLTH(Ns:Nf),'r-'); hold off axis([min(to)-0.25 max(to)+0.25... min([min(WLTH(Nf:Ns)),min(u_nm(:,4)),min(u_op(2:Nu,4))])-0.3,... max([max(WLTH(Nf:Ns)),max(u_nm(:,4)),max(u_op(2:Nu,4))])+0.3]); title('Household Wealth WLTH'); ylabel('[ billions of 1982 USD ]'); grid legend('nominal','optimal','original',-1) % subplot(3,3,6) % axis('off') % text(0.1,0.75,'L e g e n d') % text(0.1,0.60,'( solid ) historical') % text(0.1,0.50,'( dotted ) nominal'); % text(0.1,0.40,'( dashed ) simulated'); % text(0.1,0.30,'( stars ) optimal'); % hold off subplot(3,3,7) plot(to,z(2,:)','g'); hold off axis([min(to)-0.25 max(to)+0.25 min(z(2,:))-3 max(z(2,:))+3]); title('Growth Rate of Oil Price POIL'); ylabel('[ percent per year ]'); grid subplot(3,3,8) plot(to,z(3,:)','g'); hold off axis([min(to)-0.25 max(to)+0.25 min(z(3,:)')-3 max(z(3,:)')+5]); title('Net Exports XM'); ylabel('[ billions of 1982 USD ]'); grid subplot(3,3,9) plot(to,z(4,:)','g'); hold off axis([min(to)-0.25 max(to)+0.25 min(z(4,:)')-15 max(z(4,:)')+15]); title('Potential Gross National Product GNPP'); ylabel('[ billions of 1982 USD ]'); grid format disp(' '); disp(' '); disp(' Diagonal of Matrix Q'); disp(' GNP YD TAX C INR IR '); disp(diag(Q(1:6,1:6))') disp(' IIN RS RL INF WINF UR '); disp(diag(Q(7:12,7:12))') disp(' '); disp(' Diagonal of Matrix R'); disp(' G RGM TR WLTH '); disp(diag(R)') format disp(' '); disp(' '); disp(' Start and Finish Time of Simulation'); disp(' ***********************************'); disp([' t_s = ', num2str(t(Ns))]) disp([' t_f = ', num2str(t(Nf))]) disp(' ');