The way I did it was by writing my own function.
%% ARX_AR_Error: function description
function [outputs] = ARX_AR_Error(c, S, beta, rho_Reg, rho_Error)
% c is the sample length
% S regressor variable
% rho_Reg the rho of the main regression
% rho_Error the rho of the error regression
M_1 =rho_Reg*diag(ones(c-1,1),-1);
M_2=zeros(c,1);
M_2(1,1)=rho_Reg;
ey_M_1=inv(eye(c)-M_1);
% Generate the AR error
MdlY = arima('AR',{rho_Error},'Constant',0,'Variance',1);
E = simulate(MdlY,c);
outputs=ey_M_1*(S*beta+E+M_2*randn(1));
Therefore, in simulating variables I simply wrote the following code
c=100;
s=randn(c,1);
rho_Reg=0.2;
rho_Error=0.3;
beta=0.5;
Y=ARX_AR_Error(c, s, beta, rho_Reg, rho_Error)
plot(Y)
