clear;

% Define parameters
N = 2500000;

% Generate random samples from a uniform distribution
z_k = rand(1, N) - 0.5;

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

% Define initial values
P_0 = diag(100*ones(4,1));
b_0_dash = [0;0;0;0];
b_star = [4;3;2;1];

% P_k
P_k = cell(1, N);
P_k{1} = P_0;
for k = 2:N
    P_k{k} = zeros(4, 4);
end

% b_k
b_k_dash = cell(1, N);
b_k_dash{1} = b_0_dash;
for k = 2:N
    b_k_dash{k} = zeros(4, 1);
end

v_k = zeros(1,length(N));
y_k = zeros(1,length(N));

for k = 1:N
    v_k(k) = b_star(1)*u_k(k+3) + b_star(2)*u_k(k+2) + b_star(3)*u_k(k+1) + b_star(4)*u_k(k);
    y_k(k) = v_k(k) + z_k(k);
    fi_k =  [u_k(k+3); u_k(k+2); u_k(k+1); u_k(k)];
    
    if k == 1
        P_prev = P_k{k};
        b_prev = b_k_dash{k};
    else
        P_prev = P_k{k-1};
        b_prev = b_k_dash{k-1};
    end
    
    P_k{k} = P_prev - ((P_prev * fi_k * (fi_k') * P_prev)/ (1 + (fi_k') * P_prev * fi_k));
    b_k_dash{k} = b_prev + P_k{k} * fi_k * (y_k(k) - fi_k' * b_prev);
end

N_values = 1000:1000:N;
b0_values = zeros(size(N_values));
b1_values = zeros(size(N_values));
b2_values = zeros(size(N_values));
b3_values = zeros(size(N_values));

for i = 1:length(N_values)
    index = N_values(i);
    b0_values(i) = b_k_dash{index}(1);
    b1_values(i) = b_k_dash{index}(2);
    b2_values(i) = b_k_dash{index}(3);
    b3_values(i) = b_k_dash{index}(4);
end

figure('Position', [100, 100, 1200, 1000]);

subplot(2, 2, 1);
plot_subplot(N_values, b_star(1), b0_values, 'Estimation of true value b_0^*=4');

subplot(2, 2, 2);
plot_subplot(N_values, b_star(2), b1_values, 'Estimation of true value b_1^*=3');

subplot(2, 2, 3);
plot_subplot(N_values, b_star(3), b2_values, 'Estimation of true value b_2^*=2');

subplot(2, 2, 4);
plot_subplot(N_values, b_star(4), b3_values, 'Estimation of true value b_3^*=1');

function plot_subplot(N_values, true_value, estimated_values, title_str)
    plot(N_values, ones(size(N_values)) * true_value, 'LineWidth', 2, 'DisplayName', 'True Value', 'Color', '#EDB120');
    hold on;
    plot(N_values, estimated_values, 'LineWidth', 1.5, 'DisplayName', 'Estimated Value', 'Color', '#7E2F8E');
    grid on;
    legend('show', 'Location', 'best', 'FontSize', 14);
    title(title_str);
    ylabel('Values');
    xlabel('Number of samples N');
end