clear all;

% Paremeters
N = 100000;
h = 0.03125;

% Generate random input u_k
u_k = 2 * rand(1, N) - 1;

% Compute w_k based on u_k
w_k = 1 - u_k.^2;

% Initialize v_k with the first element and compute the rest
v_k = zeros(1, N);
v_k(1) = 1;
for i = 2:N
    v_k(i) = 0.5 * v_k(i-1) + w_k(i);
end

% Generate output y_k based on v_k
y_k = v_k + rand(1, N) - 0.5;

% Generate 201 numbers grid from -1 to 1
u = linspace(-1, 1, 201);
R = zeros(1,length(u));
theoretical = 1 - u.^2;
theoretical_shifted = theoretical + mean(theoretical);

for i = 1:length(u)
    nom = 0;
    den = 0;
    for k = 1:N
        nom = nom + y_k(k)*rectangular_kernel(u_k(k), u(i), h);
        den = den + rectangular_kernel(u_k(k), u(i), h);
    end
    R(i) = nom/den;
end

figure;
plot(u, theoretical_shifted, 'LineWidth', 4,'DisplayName',"Nonlinear Characteristic of Static Component (Shifted)");
hold on;
plot(u, theoretical, 'LineWidth', 4,'DisplayName',"Nonlinear Characteristic of Static Component");
hold on;
plot(u, R, 'LineWidth', 2,'DisplayName','Estimated Regression');
legend('show',Location="best");
title(['Estimated and Theoretical Values for h=', num2str(h)]);
grid on;

% Function Definition
function result = rectangular_kernel(u_k, u, h)
    valid_range = abs(u_k - u) <= h;
    result = valid_range / h;
end