dependent error distributions generated

I need to present in a graphic way four examples of dependent error distributions generated using an AR(1) model with different values of the autocorrelation coefficient.

Four further models are used to show the QR behavior with respect to the error independence assumption. In particular, starting from an autoregressive error term of order 1: e_AR(1)_[rho] → e_i = ρ * e_(i−1) + a_i

where ai ∼ N(μ = 0, σ = 1) and using the values ρ = {−0.2, +0.2, −0.5, +0.5} and i have the following models:

model1 → y_1 = 1 + 2x + e_AR(1)[ρ=−0.2];

model2 → y_2 = 1 + 2x + e_AR(1)[ρ=+0.2];

model3 → y_3 = 1 + 2x + e_AR(1)[ρ=−0.5];

model4 → y_4 = 1 + 2x + e_AR(1)[ρ=+0.5].

Below there is my work, i didn't use the function regARIMA('AR',{-0.2},'Variance',1); because i don't really know it

%% INPUT Configures
N = 1e5; % Number of samples
x = -5:0.001:5; % A row vector of equally spaced numbers
X = randn(1,size(x,2));
%% Plot Configuration
Label = 14;
Legend = 12;

%% Main Program
% Define theorical PDF of standard normal distribution

fNorm = @(x) 1/sqrt(2*pi) * exp(-x.^2/2);
% Generate vector of normally distribuited random numbers
% with mean 0 and variance 1
%X = randn(1,N);
Vett_e_AR = [ - 0.2; + 0.2; - 0.5; + 0.5 ];
rnd = normrnd(0,1,1,size(x,2));
e_AR7 = Vett_e_AR(1,1) + rnd(1,1);
e_AR8 = Vett_e_AR(2,1) + rnd(1,1);
e_AR9 = Vett_e_AR(3,1) + rnd(1,1);
e_AR10 = Vett_e_AR(4,1) + rnd(1,1);

for j = 2 : size(Vett_e_AR,1)
    R =  normrnd(0,1);   
    e_AR7(j) = Vett_e_AR(1,1) * e_AR7(j-1) + R;
    e_AR8(j) = Vett_e_AR(2,1) * e_AR8(j-1) + R;
    e_AR9(j) = Vett_e_AR(3,1) * e_AR9(j-1) + R;
    e_AR10(j) = Vett_e_AR(4,1) * e_AR10(j-1) + R;  
end

for i = 1 : size(rnd,2)
    y_7(i) = 1 + 2 * rnd(i) + e_AR7(1,end);
    y_8(i) = 1 + 2 * rnd(i) + e_AR8(1,end);
    y_9(i) = 1 + 2 * rnd(i) + e_AR9(1,end);
    y_10(i) = 1 + 2 * rnd(i) + e_AR10(1,end);
end
 fNorm = @(y_7) 1/sqrt(2*pi) * exp(-e_AR7.^2/2);
 plot(y_7, fNorm(y_7), 'r-', 'Linewidth', 2);
 grid on; axis tight; grid minor; hold on;
  fNorm_8 = @(y_8) 1/sqrt(2*pi) * exp(-y_8.^2/2);
  plot(x, fNorm_8(x), 'b-', 'Linewidth', 2);
  fNorm_9 = @(e_AR9) 1/sqrt(2*pi) * exp(-e_AR9.^2/2);
  plot(x, fNorm_9(x), 'g-', 'Linewidth', 2);
fNorm_10 = @(e_AR10) 1/sqrt(2*pi) * exp(-e_AR10.^2/2);
plot(x, fNorm_10(x), 'g-', 'Linewidth', 2);
hold off