Monte Carlo simulation problem

You can use the function, rndMVn to take draws from the multivariate normal distribution and if you only want parameter estimates, you can use the slash operator '/' to compute least squares estimates. Something like this would work for a model with a zero intercept:

beta_true = { 1.2, -0.5 };

//Covariance between 'x1' and 'x2'
cov = {   1 0.6,
        0.6 1 };

//Mean of 'x1' and 'x2'
mu = { 0, 0 };

//Number of rows per draw of 'x'
r = 100;

//Number of iterations of Monte Carlo simulation
iters = 1e5;

//Pre-allocate vector to contain
//parameter estimates for all iterations
beta_hat = zeros(iters, 2);

for i(1, iters, 1);
    //Take draws from multivariate normal distribution
    //parameterized by 'mu' and 'cov' from above
    x = rndMVn(r, mu, cov);
    
    //Create error term: err ~ N(0, 1)
    err = rndn(r, 1);
    
    //Create simulated 'y'
    y = x * beta_true + err;
    
    //Estimate beta_1 and beta_2
    //Notice transpose operator
    beta_hat[i,.] = (y / x)';
endfor;

After that you could draw a histogram of the estimated values of each estimated parameter:

//Draw histogram of beta_1
plotHistP(beta_hat[.,1], 30);

//Add histogram of beta_2 to plot
plotAddHistP(beta_hat[.,2], 30);

or calculate descriptive statistics, like this:

call dstatmt("", beta_hat);