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

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
Quartz | Level 8

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
Quartz | Level 8
Dear Rick

I will definitely follow the discussion.

Thank you so much.

SAS Innovate 2025: Call for Content

Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!

Submit your idea!

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
  • 1167 views
  • 2 likes
  • 2 in conversation