# 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 ;

//Calibrate

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

epfun=pstar;
eyfun=0.01;

// Model equationsh

model;

// 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;

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));

end;

initval;
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;
derrp=0.1;
gg=0.01;
end;

resid;

% initval; %this is for numeric steady state
% steady_state_model; %this is for analitic ss
%initializarea valorilor in ambele cazuri
%end

shocks;
var epsilon; stderr sigmay;
var eta; stderr sigmap;
var u; stderr sigmar;
end;

%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);
model_diagnostics;
%identification;
%resid;

% 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, :);
end
``````

Problem:
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?