Hierarchical Age-Period-Cohort Variance Function Regression Models Sample Codes:
                                                                                                                                   
************* step1 regression: model for the mean (Beta)-------------------------------;                                                                          
proc mixed data=temp2 covtest CL;                                                                                                       
        class year cohort;                                                                                                              
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=a;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp3;                                                                                                                             
set a;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
                                                                                                                                        
                                                                                                                                        
************* step2 regression: model for the variance (Gamma)-------------------------------;                                                                          
proc glimmix data=temp3 noitprint;                                                                                                      
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out1 pred=pred1;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp4;                                                                                                                             
set out1;                                                                                                                               
weight1=1/exp(pred1);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;                                                                                                                                    
                                                                                                                                        
                                                                                                                                        
/*iterate until the results are converged. I only list 5 iterations here. You may need                                                
more iterations to converge the results*/                                                                                               
                                                                                                                                        
/*iteration 1*/                                                                                                                         
proc mixed data=temp4 covtest CL;                                                                                                       
        class year cohort;                                                                                                              
        weight weight1;                                                                                                                 
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=b;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp5;                                                                                                                             
set b;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
proc glimmix data=temp5 noitprint;                                                                                                      
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out2 pred=pred2;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp6;                                                                                                                             
set out2;                                                                                                                               
weight2=1/exp(pred2);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;                                                                                                                                    
                                                                                                                                        
/*iteration 2*/                                                                                                                         
proc mixed data=temp6 covtest CL;                                                                                                       
        class year cohort;                                                                                                              
        weight weight2;                                                                                                                 
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=c;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp7;                                                                                                                             
set c;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
proc glimmix data=temp7 noitprint;                                                                                                      
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out3 pred=pred3;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp8;                                                                                                                             
set out3;                                                                                                                               
weight3=1/exp(pred3);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;                                                                                                                                    
                                                                                                                                        
/*iteration 3*/                                                                                                                         
proc mixed data=temp8 covtest CL;                                                                                                       
        class year cohort;                                                                                                              
        weight weight3;                                                                                                                 
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=d;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp9;                                                                                                                             
set d;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
proc glimmix data=temp9 noitprint;                                                                                                      
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out4 pred=pred4;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp10;                                                                                                                            
set out4;                                                                                                                               
weight4=1/exp(pred4);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;                                                                                                                                    
                                                                                                                                        
/*iteration 4*/                                                                                                                         
proc mixed data=temp10 covtest CL;                                                                                                      
        class year cohort;                                                                                                              
        weight weight4;                                                                                                                 
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=e;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp11;                                                                                                                            
set e;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
proc glimmix data=temp11 noitprint;                                                                                                     
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out5 pred=pred5;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp12;                                                                                                                            
set out5;                                                                                                                               
weight5=1/exp(pred5);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;                                                                                                                                    
                                                                                                                                        
/*iteration 5*/                                                                                                                         
proc mixed data=temp12 covtest CL;                                                                                                      
        class year cohort;                                                                                                              
        weight weight5;                                                                                                                 
        model y = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution outpred=f;                                                                  
        random year cohort /solution;                                                                                                   
        title "Model3: cohortort and period effects model";                                                                             
run;                                                                                                                                    
data temp13;                                                                                                                            
set f;                                                                                                                                  
res=resid*resid;                                                                                                                        
run;                                                                                                                                    
proc glimmix data=temp13 noitprint;                                                                                                     
        class year cohort;                                                                                                              
        model res = age1_c age2_c x1 x2 x3 x4 x5 x6 /solution                                                                           
          dist=gamma link=log;                                                                                                          
        random year cohort /solution;                                                                                                   
        NLOPTIONS TECHNIQUE=NMSIMP;                                                                                                     
        output out=out6 pred=pred6;                                                                                                     
        title "Model4: cohortort and period effects model of residuals";                                                                
run;                                                                                                                                    
data temp14;                                                                                                                            
set out6;                                                                                                                               
weight6=1/exp(pred6);                                                                                                                   
drop pred stderrpred df alpha lower upper resid;                                                                                        
run;