SAS/IML Matrix Maximum Likelihood Estimation - Do loop?

[LZHAO006]

Hi, I need to solve the maximum likelihood estimation problem. Here is the equation

all the matrix and scalar are given except the highlighted gamma*_hat (1x1 matrix), Beta*_hat (6x1 matrix) and Sigma*_hat (6x6 matrix).  

Here's my part of proc iml codes for the above equations. (N = 6, T =60)

I already have the initial guess number for beta and sigma. I need to plug my initial guess beta and sigma to gamma equation and get a new gamma, and then use this new gamma to calculate a new beta and sigma, and repeat this process until the results convergence. 

	/*gamma*/
		gamma = ((I - beta)`*inv(sigma)*(mu - beta*mu_m))/ ((I - beta)`*inv(sigma)*(I - beta));
	/*beta*/
			/**Numerator*/
		Nu = j(N,1,0);	D = 0;
		do i = 1 to T by 1;
			Nu = Nu + (R[i, ]` - gamma*ID)*(Rm[i,1] - gamma);
		end;
			/*Denominator*/
		do i = 1 to T by 1;
			D = D + (Rm[i,1] - gamma)##2;
		end;

		beta = Nu/D;
	
	/*sigma*/
	sigma = j(6,6,0);
		do i = 1 to T by 1;
			sigma = sigma + (Rt[i, ]`-gamma*(I-beta) - beta*Rm_b[i,1])*(Rt[i, ]`-gamma*(I-beta) - beta*Rm[i,1])`;
		end;
	sigma_b = sigma_b/T;

Does anyone know how can I repeat these codes? should I use do loop? is there a repeat code function in iml?  I really appreciate any help! 

[pink_poodle]

This article might be relevant to your question:

https://blogs.sas.com/content/iml/2013/06/19/macros-and-loops.html

[Ksharp]

Why not post it at IML forum @Rick_SAS is there .also Rick wrote several blogs about it before . you could search MLE at his blog.

[Rick_SAS]

1. What are you trying to accomplish? Are you estimating parameters in a mixed model? There might be better ways to solve this problem in SAS, so let us know what your ultimate goal is.

2. Could you post values for 

I, mu, mu_m, Rt and Rm?

Also, in your SAS code you refer to Rm_b, which doesn't seem to appear in the equations.

[LZHAO006]

Here's the post value of  Rt

[LZHAO006]

HI Rick, my ultimate goal is to estimate the gamma*, beta* and sigma* using the provided equation. Right now I already have the initial estimate value for beta and sigma. and I would like to run a loop (maybe 1000times?) for my equation codes to get the gamma*, beta* and sigma*.

This is part of the Black version of CAPM calculation. so I add all value name with "_b" to remind me this is the black version in my original code. Sorry I didn't delete when I posted my questions.

Please see the attachments of my post values and initial guest for beta and sigma. I is the 6x1 identity matrix. I added in my code. Thank you very much!

I = j(6,1,1);
T = 60;

/*gamma*/
		gamma = ((I - beta)`*inv(sigma)*(mu - beta*mu_m))/ ((I - beta)`*inv(sigma)*(I - beta));

/*beta*/
			/**Numerator*/
		Nu = j(N,1,0);	D = 0;
		do i = 1 to T by 1;
			Nu = Nu + (Rt[i, ]` - gamma*ID)*(Rm[i,1] - gamma);
		end;
			/*Denominator*/
		do i = 1 to T by 1;
			D = D + (Rm[i,1] - gamma)##2;
		end;

		beta = Nu/D;
	
/*sigma*/
	sigma = j(6,6,0);
		do i = 1 to T by 1;
			sigma = sigma + (Rt[i, ]`-gamma*(I-beta) - beta*Rm[i,1])*(Rt[i, ]`-gamma*(I-beta) - beta*Rm[i,1])`;
		end;
	sigma = sigma/T;

[Rick_SAS · Accepted solution]

SAS is not a case-sensitive language. Thus your matrix I and the iteration variable i are the same variable. As soon as you execute the first loop, the I matrix is overwritten/replaced.

Do know from theory that the iterative computation of those equations will converge to a solution? If so,

write a loop such as the following (untested, but hopefully you can see the main idea):

beta = GUESS1; gamma = GUESS2; SIGMA = GUESS3;
convergence = 0; /* false */
do k = 1 to 500 until (convergence);
  beta_old = beta; gamma_old = gamma; Sigma_old = Sigma;
   /* compute/update beta, gamma, and Sigma according to formulas */
  beta_diff = norm(abs(beta - beta_old));   /* or use max(abs( ... )) */
  gamma_diff = norm(abs(gamma - gamma_old));
  Sigma_diff =  norm(abs(Sigma - Sigma_old));
  if beta_diff <= 1e-3 & gamma_diff <= 1e-3 & Sigma_diff <= 1e-3 then
      convergence=1;
  if mod(k, 25)=0 then
      print k beta_diff gamma_diff Sigma_diff;
end;

print converged k beta_diff gamma_diff Sigma_diff;

[LZHAO006]

It works! thank you very much!!!