BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Lyson
Obsidian | Level 7

Helo family!

 

I am looking for a code to specify arh(1) covariance structure in the PROC IML simulation procedure. I saw the IML code for specifying the R-side ar(1)  covariance structure by Gibbs & Kiernan (2020), but I am struggling to figure out how to capture arh(1) instead and generate my simulation samples. Here is the code below by Gibbs & Kiernan (2020), when the R-side covariance structure is ar(1):

proc iml worksize=100;
  resid=2;
  rho=0.8;
  m=5; /* m = number of time points  */
  cov=j(m,m,1);
  do i=1 to m;
      do j=1 to m;
        cov=[i, j]=resid*rho**abs(i-j);
    end;
end; mean={0 0 0}; n=200; call randseed(61345); z=RandNormal(n, mean, cov); create AR1Errors from z; append from z; quit;

Kindly assist with the parameters (code) required for arh(1) structure in a similar code as given above. Alternative approaches (codes) are most welcome.

 

Thank you for your assistance.

 

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

The ARH(1) structure is shown in the PROC MIXED documentation.

The structure is the Hadamard product of an outer-product matrix sigma*sigma` and an AR1 matrix (which is a Toeplitz matrix). I think the following IML program will sample normal errors from an ARH1 covariance matrix:

 

proc iml;
sigma = {8, 16, 4, 8, 1};
rho=0.5;
dim = nrow(sigma);

/* AR1 correlation structure:
   https://blogs.sas.com/content/iml/2012/11/05/constructing-common-covariance-structures.html
*/
start AR1Corr(dim, rho);
   u = cuprod(j(1,dim-1,rho)); /* cumulative product */
   return( toeplitz(1 || u) );
finish;

/* ARH(1) is Hadamard product of outer product sigma*sigma` 
   and the AR(1) correlation matrix */
OutProd = sigma * sigma`;
AR1 = AR1Corr(dim, rho);
ARH1Cov = OutProd # AR1;
print OutProd, AR1, ARH1Cov;

mean = j(1, dim, 0);
n = 200;
call randseed(61345);
z = RandNormal(n, mean, ARH1Cov);

create ARH1Errors from z[c=('x1':'x5')];
append from z;
quit;

/* visualize */
proc sgscatter data=ARH1Errors;
matrix x1-x5;
run;

View solution in original post

4 REPLIES 4
Rick_SAS
SAS Super FREQ

The ARH(1) structure is shown in the PROC MIXED documentation.

The structure is the Hadamard product of an outer-product matrix sigma*sigma` and an AR1 matrix (which is a Toeplitz matrix). I think the following IML program will sample normal errors from an ARH1 covariance matrix:

 

proc iml;
sigma = {8, 16, 4, 8, 1};
rho=0.5;
dim = nrow(sigma);

/* AR1 correlation structure:
   https://blogs.sas.com/content/iml/2012/11/05/constructing-common-covariance-structures.html
*/
start AR1Corr(dim, rho);
   u = cuprod(j(1,dim-1,rho)); /* cumulative product */
   return( toeplitz(1 || u) );
finish;

/* ARH(1) is Hadamard product of outer product sigma*sigma` 
   and the AR(1) correlation matrix */
OutProd = sigma * sigma`;
AR1 = AR1Corr(dim, rho);
ARH1Cov = OutProd # AR1;
print OutProd, AR1, ARH1Cov;

mean = j(1, dim, 0);
n = 200;
call randseed(61345);
z = RandNormal(n, mean, ARH1Cov);

create ARH1Errors from z[c=('x1':'x5')];
append from z;
quit;

/* visualize */
proc sgscatter data=ARH1Errors;
matrix x1-x5;
run;
Lyson
Obsidian | Level 7

Hi Rick,

 

Thank you so much. You are very right. I got a chance to go through the suggested code and additional reading material. It is quite clear. I thank you.

Lyson
Obsidian | Level 7
Dear Rick

I will definitely follow the discussion.

Thank you so much.

Ready to join fellow brilliant minds for the SAS Hackathon?

Build your skills. Make connections. Enjoy creative freedom. Maybe change the world. Registration is now open through August 30th. Visit the SAS Hackathon homepage.

Register today!
Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 4 replies
  • 900 views
  • 2 likes
  • 2 in conversation