NaN in squaring a difference, basic NK model

I have the following dynare code/model:

// Endogenous variables
var y p r ep_extrap eps ey_extrap eys CRp FRp CRy FRy alfap alfay derrp gg;
// Exogenous variablesee
varexo epsilon eta u;
parameters mm pstar a1 a2 b1 b2 c1 c2 c3 sigmay sigmap sigmar epfun eyfun rho rhoBH ;
mm = 2;         %switching parameter gamma in Brock Hommes % ? din ec 12 si 13 (n-ai nevoie sa-l modifici)
pstar = 0.01;      % the central bank's inflation target
a1 = 0.5;       %coefficient of expected output in output equation
a2 = -0.4;      %a is the interest elasticity of output demand
b1 = 0.5;       %b1 is coefficient of expected inflation in inflation equation
b2 = 0.05;      %b2 is coefficient of output in inflation equation
c1 = 1.5;       %c1 is coefficient of inflation in Taylor equation
c2 = 0.5;       %c2 is coefficient of output in Taylor equation
c3 = 0.5;       %interest smoothing parameter in Taylor equation
sigmay = 0.1;   %standard deviation shocks output 
sigmap = 0.1;   %standard deviation shocks inflation
sigmar = 0.1;   %standard deviation shocks Taylor
mu=1;           %adaptive paramneter in fundamnetalist rule
rho=0.5;        %rho in mean squares errors ec 9' de la nota de subsol
rhoBH=0.0;      %rho in Brock-Hommes switching 
// Model equationsh
    // NK equations
    y = a1*eys + (1-a1)*y(-1) + a2*(r-eps) + epsilon;
    p = b1*eps + (1-b1)*p(-1) + b2*y + eta;
    r= c1*p + c2*y + c3*r(-1) + u;  
    // Expectations updates
    eps = alfap*ep_extrap + (1-alfap)* epfun;  
    eys = alfay*ey_extrap + (1-alfay)*eyfun; 
    ep_extrap = p(-1); 
    ey_extrap = y(-1);      %AR1 extrapolators
    derrp = ep_extrap - p;
    gg =  derrp*derrp;
    // Cost ratios and preference updates
    CRp = rho * CRp(-1) - (1-rho) * (ep_extrap - p)^2;
    FRp = rho * FRp(-1) - (1-rho) * (epfun - p)^2;
    CRy = rho * CRy(-1) - (1-rho) * (ey_extrap - y)^2;
    FRy = rho * FRy(-1) - (1-rho) * (eyfun - y)^2;
    alfap = rhoBH*alfap(-1) + (1-rhoBH)*exp(mm*CRp)/(exp(mm * CRp) + exp(mm * FRp)); 
    alfay = rhoBH*alfay(-1) + (1-rhoBH)*exp(mm*CRy)/(exp(mm * CRy) + exp(mm * FRy));
y = 0;
p = 0.01;
r = 0;
ep_extrap = 0;
ey_extrap = 0.1;
eps = 0;
eys = 0;
CRp = 0;
FRp = 5;
CRy = 10;
FRy = 6;
alfap = 0.5;
alfay = 0.5;
% initval; %this is for numeric steady state
% steady_state_model; %this is for analitic ss
%initializarea valorilor in ambele cazuri
var epsilon; stderr sigmay;
var eta; stderr sigmap;
var u; stderr sigmar;
%varobs y p r;
// Check model specification 
stoch_simul(order=1, irf=20, periods=200, solve_algo=4);
%stoch_simul(order=1, irf=20, periods=200, drop=10);
% global M_ oo_; % Make sure global variables are accessible
% indices_derrp = strmatch('derrp', M_.endo_names, 'exact'); % Get the index for derrp
% indices_gg = strmatch('gg', M_.endo_names, 'exact'); % Get the index for gg
% oo_.endo_simul(indices_gg, :) = oo_.endo_simul(indices_derrp, :).^2; % Update gg
% Initialize a struct to hold the simulated data
simulated_data = struct();
% Loop through all endogenous variable names
for i = 1:length(M_.endo_names)
    % Get the variable name
    var_name = M_.endo_names{i};
    % Extract the simulated data for this variable
    simulated_data.(var_name) = oo_.endo_simul(i, :);

I obtain NaNs for CRp and CRy.

Let’s take only CRp as an example.
The problem comes from squaring the difference between ep_extrap and p:
(ep_extrap - p)^2

I calculated separately the difference and it seems OK, but when dynare squares it, it produces NaNs. The exact same thing happens for CRy.

The simulations are OK and I can’t find any extreme values (either positive or negative) to lead to this.

Does anyone have any idea on how to get rid of the issue of obtaining NaNs?
Is there anything wrong in the declaration of my model?

Thanks in advance!

In Dynare 6 I obtain zeros for those variables, which is consistent with the depicted policy rules that show no shock affecting these variables.

Did you run the code I provided and you obtain zeros instead of NaNs?

That’s weird, I also use Dynare 6 and I get NaNs.

Did you do any changes to the code or would you recommend me modifying the existing code in some way?


I obtain NaN for variance decompositions. But again, that is expected if the variance is 0.