How to nest temporary shock model(rho=0.01) and permanent shock model(rho=0.99) in one graph for rest similar data?

I am using following code even by changing the model_type on top i am not able to get it in one graph.

%%%%%%%%%%%
% A Simplified Template for solving linear rational expectation models
%%%%%%%%%%%
% example: the neoclassical growth model with complete depreciation
model_type =1; %choose the model type (1 or 2) 
%% Calibrate the model
alpha = 0.33;
beta = 0.99;
delta = 0.025;
if model_type==2
    rho=0.01;
else 
    model_type = 1;
    rho=0.99;                  %technology shock
end
gamma = 1;
gamma_N = 1;
ssa = 1;                %steady state technology
ssn = 1;                %steady state labor
%% Calculate the initial steady state values
k = ssn*((1-beta*(1-delta))/(beta*alpha*ssa))^(1/(alpha-1));       %S.S euler equation
i = delta*k;                                                       %S.S law of motion
y = k^(alpha);                                                     %S.S production function
c = y-i;                                                           %S.S ARC equation
r = alpha*ssa*(k)^(alpha-1)*(ssn)^(1-alpha);                       %S.S labor demand equation
w = (1-alpha)*ssa*k^(alpha)*ssn^(-alpha);                          %S.S capital demand equation
theta = c^(1-alpha)*ssa*(k)^(alpha)*ssn^(-alpha-gamma_N);          %S.S labor supply equation
c_y = c/y;                                                         %consumption output ratio
k_y = k/y;                                                         %capital output ratio
%% Define variable numbers
nc = 1;
nk = 2;
nn = 3;
ny = 4;
nw = 5;
nr = 6;
ni = 7;
na = 8;
% Define the matrices and relevant paramters
nvar = 8;          %number of variables
nshocks =1;        %number of shocks
nfore = 2;         %number of forecast errors
g0 = zeros(nvar,nvar);
g1 = zeros(nvar,nvar);
psi = zeros(nvar,nshocks);
pi = zeros(nvar,nfore);
cc = zeros(nvar,1);
%% Input the equations
%%%%%%%%%
% (1) Euler Equation
%%%%%%%%%
g0(1,nc) = -gamma;
g0(1,nr) = beta*r;
g1(1, nc) = -gamma;
pi(1,1) = -gamma;
pi(1,2) = beta*r;
%%%%%%%%%%%
% (2) Law of motion of caliptal
%%%%%%%%%%%
g0(2,nk) = 1;
g0(2, ni) = -delta;
g1(2,nk) = 1-delta;
%%%%%%%%%%%%
% (3) labor supply equation
%%%%%%%%%%%%
g0(3,nn) = gamma_N;
g0(3, nc) = gamma;
g0(3,nw) = -1;
%%%%%%%%%%%%
% (4) Production function
%%%%%%%%%%%%
g0(4,nn) = -(1-alpha);
g0(4,ny) = 1;
g0(4,na) = -1;
g1(4,nk) = alpha;
%%%%%%%%%%%%
% (5) labor demand equation
%%%%%%%%%%%%
g0(5,nn) = alpha;
g0(5, nw) = 1;
g0(5,na) = -1;
g1(5,nk) = alpha;
%%%%%%%%%%%%
% (6) capital demand equation
%%%%%%%%%%%%
g0(6,nn) = -(1-alpha);
g0(6, nr) = 1;
g0(6,na) = -1;
g1(6,nk) = (alpha-1);
%%%%%%%%%%%%
% (7) ARC
%%%%%%%%%%%%
g0(7,nc) = c;
g0(7, ni) = i;
g0(7,ny) = -y;
%%%%%%%%%%%%
% (8) stochastic process for technology
%%%%%%%%%%%%
if model_type==2
    g0(8,na) =1;
    psi(8,1) =1;
else
    model_type = 1;
    g0(8,na) =1;
    g1(8,na) =rho;
    psi(8,1) =1;
end
%% solve the model
[G,C0,M,fmat,fwt,ywt,gev,eu]=gensys(g0,g1,cc,psi,pi);

Unrecognized function or variable 'gensys'.

% Check if existence and uniqueness
eu;
if eu(1,1) ==1 && eu(2,1)==1
    disp([' eu=    ' mat2str(eu)])
    disp('There exists a unique equilibrium.')
else
    disp([' eu=    ' mat2str(eu)])
    return
end
%% Generate some impulse response functions
% Size of the Shock to each Random Variable
sigma_e = 1;   % 1% Shock to e
% --- Initial Response  ---------------------------------------------------
periods=60; %Length of the Response
IMF=zeros(nvar,nshocks,periods);  
sigmamatrix=sigma_e;
IMF(:,:,1)= M*sigmamatrix;  
% --- Response Assuming no Additional Shocks  in the Future ---------------
for s=2:periods;
    IMF(:,:,s)= G * IMF(:,:,s-1);
end
% --- Eliminate Very Small Responses --------------------------------------
for i = 1:nvar
    for j = 1:nshocks
        for k = 1:periods
            if abs(IMF(i,j,k))<10^-10
                IMF(i,j,k) = 0;
            end
        end
    end
end
% --- Plot Response of Each Variable --------------------------------------
Years = (0:(periods-1))/4;  % converts time dimension to years
Quarters = 0:periods-1;
varnames = {'Consumption';'Capital';'labor';'output';'wage';'interest';'investment';'technology'}; 
set(0,'DefaultTextFontsize',20, ...
    'DefaultTextFontWeight','normal', ...
    'DefaultAxesFontsize',14)
for k = 1:nshocks
    figure(k);
    for j = 1:nvar
        subplot(3,3,j);
        IMF1(1,:)=IMF(j,k,1:periods);
        plot(Quarters,IMF1,'b-','LineWidth',2,'MarkerSize',3);
        grid on;
        title(char(varnames(j)), 'fontsize', 14);
    end
end
return

[Avni Agrawal]

Hi, i'm trying to run the code and getting error at line 109: "Unrecognized function or variable 'gensys'."

Can you share this function code as well or is this function the part of any toolbox?

Answer 1

[Walter Roberson]

http://sims.princeton.edu/yftp/gensys/mfiles/

as described at http://sims.princeton.edu/yftp/Macro99/UsingGensys.pdf

Answer 2

[Avni Agrawal]

Hi,

I understand that you aim to create sub-plots for both model types. To accomplish this, you can segment the model logic from the calculation of initial steady state values to the elimination of very small responses.

Your function might resemble the following:

function [periods, nshocks, nvar, IMF] =  shock_model(alpha, beta, delta, ssa, ssn,gamma, gamma_N, rho, model_type)
    %% Calculate the initial steady state values
    k = ssn*((1-beta*(1-delta))/(beta*alpha*ssa))^(1/(alpha-1));       %S.S euler equation
    i = delta*k;                                                       %S.S law of motion
    y = k^(alpha);                                                     %S.S production function
    c = y-i;                                                           %S.S ARC equation
    r = alpha*ssa*(k)^(alpha-1)*(ssn)^(1-alpha);                       %S.S labor demand equation
    w = (1-alpha)*ssa*k^(alpha)*ssn^(-alpha);                          %S.S capital demand equation
    theta = c^(1-alpha)*ssa*(k)^(alpha)*ssn^(-alpha-gamma_N);          %S.S labor supply equation
    c_y = c/y;                                                         %consumption output ratio
    k_y = k/y;                                                         %capital output ratio
    %% Define variable numbers
    nc = 1;
    nk = 2;
    nn = 3;
    ny = 4;
    nw = 5;
    nr = 6;
    ni = 7;
    na = 8;
    % Define the matrices and relevant paramters
    nvar = 8;          %number of variables
    nshocks =1;        %number of shocks
    nfore = 2;         %number of forecast errors
    g0 = zeros(nvar,nvar);
    g1 = zeros(nvar,nvar);
    psi = zeros(nvar,nshocks);
    pi = zeros(nvar,nfore);
    cc = zeros(nvar,1);
    %% Input the equations
    %%%%%%%%%
    % (1) Euler Equation
    %%%%%%%%%
    g0(1,nc) = -gamma;
    g0(1,nr) = beta*r;
    g1(1, nc) = -gamma;
    pi(1,1) = -gamma;
    pi(1,2) = beta*r;
    %%%%%%%%%%%
    % (2) Law of motion of caliptal
    %%%%%%%%%%%
    g0(2,nk) = 1;
    g0(2, ni) = -delta;
    g1(2,nk) = 1-delta;
    %%%%%%%%%%%%
    % (3) labor supply equation
    %%%%%%%%%%%%
    g0(3,nn) = gamma_N;
    g0(3, nc) = gamma;
    g0(3,nw) = -1;
    %%%%%%%%%%%%
    % (4) Production function
    %%%%%%%%%%%%
    g0(4,nn) = -(1-alpha);
    g0(4,ny) = 1;
    g0(4,na) = -1;
    g1(4,nk) = alpha;
    %%%%%%%%%%%%
    % (5) labor demand equation
    %%%%%%%%%%%%
    g0(5,nn) = alpha;
    g0(5, nw) = 1;
    g0(5,na) = -1;
    g1(5,nk) = alpha;
    %%%%%%%%%%%%
    % (6) capital demand equation
    %%%%%%%%%%%%
    g0(6,nn) = -(1-alpha);
    g0(6, nr) = 1;
    g0(6,na) = -1;
    g1(6,nk) = (alpha-1);
    %%%%%%%%%%%%
    % (7) ARC
    %%%%%%%%%%%%
    g0(7,nc) = c;
    g0(7, ni) = i;
    g0(7,ny) = -y;
    %%%%%%%%%%%%
    % (8) stochastic process for technology
    %%%%%%%%%%%%
    if model_type==2
        g0(8,na) =1;
        psi(8,1) =1;
    else
        model_type = 1;
        g0(8,na) =1;
        g1(8,na) =rho;
        psi(8,1) =1;
    end
    %% solve the model
    [G,C0,M,fmat,fwt,ywt,gev,eu]=gensys(g0,g1,cc,psi,pi);
    % Check if existence and uniqueness
    eu;
    if eu(1,1) ==1 && eu(2,1)==1
        disp([' eu=    ' mat2str(eu)])
        disp('There exists a unique equilibrium.')
    else
    disp([' eu=    ' mat2str(eu)])
        return
    end
    %% Generate some impulse response functions
    % Size of the Shock to each Random Variable
    sigma_e = 1;   % 1% Shock to e
    % --- Initial Response  ---------------------------------------------------
    periods=60; %Length of the Response
    IMF=zeros(nvar,nshocks,periods);  
    sigmamatrix=sigma_e;
    IMF(:,:,1)= M*sigmamatrix;  
    % --- Response Assuming no Additional Shocks  in the Future ---------------
    for s=2:periods
       IMF(:,:,s)= G * IMF(:,:,s-1);
    end
    % --- Eliminate Very Small Responses --------------------------------------
    for i = 1:nvar
       for j = 1:nshocks
          for k = 1:periods
             if abs(IMF(i,j,k))<10^-10
                IMF(i,j,k) = 0;
             end
          end
       end
    end
end

To plot the responses for both model types, you should call the function twice, passing appropriate variables. Here is the code:

% example: the neoclassical growth model with complete depreciation
model_type =2; %choose the model type (1 or 2) 
%% Calibrate the model
alpha = 0.33;
beta = 0.99;
delta = 0.025;
if model_type==2
    rho=0.01;
else 
    model_type = 1;
    rho=0.99;                  %technology shock
end
gamma = 1;
gamma_N = 1;
ssa = 1;                %steady state technology
ssn = 1;                %steady state labor

% --- Plot Response of Each Variable --------------------------------------

%call our custom defined shock_model
[periods, nshocks, nvar, IMF] = shock_model(alpha, beta,delta, ssa, ssn,gamma,gamma_N, 0.01, 2); %model type 1
[periods1, nshocks1, nvar1, IMF1] = shock_model(alpha, beta,delta, ssa, ssn,gamma,gamma_N, 0.99, 1); %model type 2

Years = (0:(periods-1))/4;  % converts time dimension to years
Quarters = 0:periods-1;
varnames = {'Consumption';'Capital';'labor';'output';'wage';'interest';'investment';'technology'}; 
set(0,'DefaultTextFontsize',20, ...
'DefaultTextFontWeight','normal', ...
'DefaultAxesFontsize',14)
for k = 1:nshocks
    figure(k);
    for j = 1:nvar
        subplot(3,3,j);
        IMF_(1,:)=IMF(j,k,1:periods);
        plot(Quarters,IMF_,'b-','LineWidth',2,'MarkerSize',3);
        hold on;
        IMF1_(1,:)=IMF1(j,k,1:periods1);
        plot(Quarters,IMF1_,'b-','LineWidth',2,'MarkerSize',3, 'Color','k');
legend('model type 1', 'model type 2');
        hold off;
        grid on;
        title(char(varnames(j)), 'fontsize', 14);
    end
end
return;

By following this approach, you can nest plots for both model types.

I hope this helps.