periods = 100;

% Parameters
sigma_epsilon = 1;
theta = 0.8;
rho_M = 0.8;
M_bar = 1;
sigma_zeta = 2;

% Initial values 
u = 0;
M = 1;

M_list = zeros(periods, 1);
y_list = zeros(periods - 1, 1);
u_list = zeros(periods, 1);


% Initialisation
M_list(1,1) = M;
u_list(1,1) = u;

% Shocks
epsilon = sigma_epsilon * randn(periods,1);
nu_M = randn(periods,1);
zeta = randn(periods,1);


for i = 2: periods
    M_list(i,1) = rho_M * M_list(i-1,1) + (1 - rho_M) * M_bar + nu_M(i,1);
    u_list(i,1) = theta * u_list(i-1,1) + sigma_zeta * zeta(i,1) + M_list(i,1) * epsilon(i,1);
end

for i = 1: periods-1
    y_list(i,1) = u_list(i+1,1) + epsilon(i+1,1); 
end

% Rename y_list 
y = y_list;
save('data_simple', 'y');

dynare basic_filter4 noclearall

% Plot different graphs showing the difference between the updated values
% and the real values
u_graph = u_list(2:100,1);
figure
plot(oo_.SmoothedVariables.u);
hold on
plot(u_graph,'r')
title('u');

M_graph = M_list(2:100,1);
figure
plot(oo_.SmoothedVariables.M);
hold on
plot(M_graph,'r')
title('M');

epsilon_graph = epsilon(2:100,1);
figure
plot(oo_.SmoothedShocks.epsilon);
hold on
plot(epsilon_graph,'r')
title('epsilon');

zeta_graph = zeta(2:100,1);
figure
plot(oo_.SmoothedShocks.zeta);
hold on
plot(zeta_graph,'r')
title('zeta');