Symmetric capital stocks

I have used a standard deterministic Solow model (example of Johannes Pfeifer) and included (in addition to the physical capital stock K) a second capital stock H in a symmetric way: CES(L,CES(K,H)). The share of both capital stocks are equal and investments are also fixed at 50%-50%.

However, if I shock the capital stocks separately by the same amount, the solutions are different, although they should be symmetric.

Another troubling issue is that the capital stock H is equal to 0 in the initial period if I shock it.

I can easily delete the endval-values except for the exogenous variables and get the same solutions (as described by Johannes Pfeifer in http://www.dynare.org/phpBB3/viewtopic.php?f=1&t=5345).

I’ve tried to find the asymmetry in the model but was not successful, and I cannot explain the asymmetric solutions to the two shocks.

Many thanks for your help!

//****************************************************************************
//Define variables
//****************************************************************************
var c           ${c}$ (long_name='consumption (intensive form)')
    k           ${k}$ (long_name='capital (intensive form)')
    h           ${h}$ (long_name='human capital (intensive form)')
    y           ${y}$ (long_name='output (intensive form)')
    invest      ${i}$ (long_name='total investment capital (intensive form)')    
    invest_h    ${i_h}$ (long_name='investment in human capital (intensive form)')
    ;
//****************************************************************************
//Define predetermined variables
//****************************************************************************
predetermined_variables k;
//****************************************************************************
//Define exogenous variables (delta_k and epsilon are time-varying in our experiment)
//****************************************************************************
varexo delta_k    ${\delta_k}$ (long_name='depreciation rate of physical capital')
       delta_h    ${\delta_k}$ (long_name='depreciation rate of human capital')
       epsilon       ${\epsilon}$ (long_name='elasticity of subsitution between K and H')
    ;
//****************************************************************************
//Define parameters
//****************************************************************************
parameters s    ${s}$ (long_name='saving rate')
    alpha       ${\alpha}$ (long_name='capital share production function')
    alpha_k     ${\alpha_k}$ (long_name='capital share in inner production function')
    s_i         ${\alpha_i}$ (long_name='share of investments in physical capital')
    delta_k_initial       ${\delta_k_0}$ (long_name='depreciation rate of physical capital')
    delta_h_initial     ${\delta_h_0}$ (long_name='depreciation rate of human capital')
    n           ${n}$ (long_name='population growth rate')
    g           ${g}$ (long_name='technology growth rate')
    epsilon_initial   ${\epsilon_0}$ (long_name='initial elasticity of substitution')
    epsilon_k    ${\epsilon_k}$ (long_name='elasticity of substitution 2')
    ;
//****************************************************************************
//set parameters
//****************************************************************************
s=0.2;
alpha=0.3;           
alpha_k=0.5;
s_i=0.5;
delta_k_initial=0.1;
delta_h_initial=0.1;
n=0;
g=0.02;
epsilon_initial=1.0001;
epsilon_k=1.0001;

        
//****************************************************************************
//enter the model equations (model-block)
//****************************************************************************
model;
[name='Law of motion capital']
(1+n+g+n * g) * k(+1)=(1-delta_k) * k+s_i * invest;
[name='Law of motion human capital']
(1+n+g+n * g) * h(+1)=(1-delta_h) * h+invest_h;
[name='resource constraint']
invest+c=y;
[name='total investments']
invest_h=(1-s_i) * invest;
[name='behavioral rule savings']
c=(1-s) * y;
[name='production function']
y=((1-alpha)+alpha * (alpha_k * k^((epsilon_k-1)/epsilon_k) + (1-alpha_k) * h^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
end;




%%%%%%%%%%%%%% Shock#2: Increase depreciation rate of physical capital %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




//****************************************************************************
//initval-block: set initial condition to steady state value with delta_k=delta_k_initial
//and delta_h=delta_h_initial        
//****************************************************************************
initval;
    epsilon=epsilon_initial;
    delta_k=delta_k_initial;
    delta_h=delta_h_initial;
    k=((1/(1-alpha))*(((delta_k+n+g+n * g)/(s_i * s))^((epsilon-1)/epsilon)-alpha * (alpha_k+(1-alpha_k) * (((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i) / (s_i)))^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1))*((epsilon-1)/epsilon))))^(epsilon/(1-epsilon));
    h=((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i) / (s_i)) * k;
    y=((1-alpha)+alpha*(alpha_k * k^((epsilon_k-1)/epsilon_k) + (1-alpha_k) * h^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
    c=(1-s) * y;
    invest=y-c;
    invest_h=(1-s_i) * invest;
end;
//****************************************************************************
//endval-block: set 
//      i)  delta_k=0.3 and epsilon=epsilon_initial and delta_h=delta_h_initial for all periods after the initial one
//      ii) terminal condition to steady state value with delta_k=0.3 and epsilon=epsilon_initial and delta_h=delta_h_initial;
//****************************************************************************
endval;
    epsilon=epsilon_initial;
    delta_h=delta_h_initial;
    delta_k=0.3;
//    k=((1/(1-alpha)) * (((delta_k+n+g+n * g)/(s_i * s))^((epsilon-1)/epsilon)-alpha * (alpha_k+(1-alpha_k) * (((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i) / (s_i))) ^ ((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon))))^(epsilon/(1-epsilon));
//    h=((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i) / (s_i)) * k;
//    y=((1-alpha)+alpha * (alpha_k * k^((epsilon_k-1)/epsilon_k) + (1-alpha_k) * h^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
//    c=(1-s) * y;
//    invest=y-c;
//    invest_h=(1-s_i) * invest;
end;
//****************************************************************************
//steady-command: compute the steady state conditional on the value of the exogenous 
//  states and use them for endval (redundant here, because we provided the steady states
//  analytically) 
//****************************************************************************
steady;
//****************************************************************************
//perfect_foresight_setup: set up the simulation for 100 periods
//you can check the settings in oo_.endo_simul (and oo_.exo_simul if there 
//were exogenous variables)
//****************************************************************************
perfect_foresight_setup(periods=100);
//****************************************************************************
//perfect_foresight_solver: compute the solution
//****************************************************************************
perfect_foresight_solver;
//****************************************************************************
//rplot command: display simulation results
//****************************************************************************

dsample 20;

rplot k h y c invest invest_h;











%%%%%%%%%%%%%% Shock#3: Increase depreciation rate of human capital %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




//****************************************************************************
//initval-block: set initial condition to steady state value with delta_h=delta_h_initial
//and delta_k=delta_k_initial            
//****************************************************************************
initval;
    epsilon=epsilon_initial;
    delta_k=delta_k_initial;
    delta_h=delta_h_initial;
    k=((1/(1-alpha)) * (((delta_k+n+g+n * g)/(s_i * s))^((epsilon-1)/epsilon)-alpha * (alpha_k+(1-alpha_k) * (((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i)/(s_i)))^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon))))^(epsilon/(1-epsilon));
    h=((delta_k+n+g+n * g)/(delta_h+n+g+n * g))*((1-s_i)/(s_i)) * k;
    y=((1-alpha)+alpha*(alpha_k * k^((epsilon_k-1)/epsilon_k) + (1-alpha_k) * h^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
    c=(1-s) * y;
    invest=y-c;
    invest_h=(1-s_i) * invest;
end;
//****************************************************************************
//endval-block: set 
//      i)  delta_h=0.3 and epsilon=epsilon_initial and delta_k=delta_k_initial for all periods after the initial one
//      ii) terminal condition to steady state value with delta_h=0.3 and epsilon=epsilon_initial and delta_k=delta_k_initial;
//****************************************************************************
endval;
    epsilon=epsilon_initial;
    delta_k=delta_k_initial;
    delta_h=0.3;
//    k=((1/(1-alpha)) * (((delta_k+n+g+n * g)/(s_i * s))^((epsilon-1)/epsilon)-alpha * (alpha_k+(1-alpha_k) * (((delta_k+n+g+n * g)/(delta_h+n+g+n * g)) * ((1-s_i) / (s_i)))^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1)) * ((epsilon-1)/epsilon))))^(epsilon/(1-epsilon));
//    h=((delta_k+n+g+n * g)/(delta_h+n+g+n * g))*((1-s_i) / (s_i)) * k;
//    y=((1-alpha)+alpha*(alpha_k * k^((epsilon_k-1)/epsilon_k) + (1-alpha_k) * h^((epsilon_k-1)/epsilon_k))^((epsilon_k/(epsilon_k-1))*((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
//    c=(1-s) * y;
//    invest=y-c;
//    invest_h=(1-s_i) * invest;
end;
//****************************************************************************
//steady-command: compute the steady state conditional on the value of the exogenous 
//  states and use them for endval (redundant here, because we provided the steady states
//  analytically) 
//****************************************************************************
steady;
//****************************************************************************
//perfect_foresight_setup: set up the simulation for 100 periods
//you can check the settings in oo_.endo_simul (and oo_.exo_simul if there 
//were exogenous variables)
//****************************************************************************
perfect_foresight_setup(periods=100);
//****************************************************************************
//perfect_foresight_solver: compute the solution
//****************************************************************************
perfect_foresight_solver;
//****************************************************************************
//rplot command: display simulation results
//****************************************************************************

dsample 20;

rplot k h y c invest ;

Because you have a timing error. It is

predetermined_variables k;

but should be

predetermined_variables k h;

Many thanks for your quick reply, it is very appreciated!