/* Program creates a vector of "price" series that are AR1.  
In the program, "numpri" is the number of individual price series.
The length of each price series is "sampsiz".  
The "numbri" AR1 coefficients are drawn from a uniform distribution,
with mean of "meanr" and standard deviation of "std".
The series have normally distributed shocks, with a standard
deviation of "sterr".
The program averages the price series, then estimates an AR1 on the average.
This is repeated "rep" times. 
The comments below indicate the assumptions on each of these parameters
in the simulations used by Imbs et. al.

Written by Charles Engel, August 2003 
Modified by Shiu-Sheng Chen, February 2004 */

rep = 2500;

numpri = 100;  /* Number of individual price series.  Imbs et. al.
                  average over 20 goods (and run panels with 10 average
                  prices.  This program only does univariate estimation,
                  calculating AR1 averaged over numpri prices.) */

sampsiz = 500;  /* Length of each time series. Imbs et. al. assume 100.*/

meanr = .94;      /* Average AR1 coefficient.  
                      Imbs et. al. take a number around .94.*/

rng = .059;         /* The range of the AR1 coefficients is plus or minus
                       rng. I believe Imbs et. al. assume rng is .15.
                       They say the variance is .0075.  For a uniform,
                       if you let x be the deviation from the mean,
                       then the variance is given by (1/3)*(range^2),
                       where the range is the maximum minus the minimum.
                       The solution to (1/12)*(range^2)=.0075
                       is range = .3, so rng = .15. */



sterr = 1;        /* The standard error of innovations to prices.  Imbs
                     et. al. set this equal to 1. */

coeffs_a = 0;  /* This vector will store the AR1 coefficients for the average
                   of prices */  
coeffs_b = 0;

meanrho = 0;   /*  This vector will store the average AR1 coefficient from 
                  the "numbri" price series -- this should equal "meanr".*/

j=1;

do until j == rep+1;

j;

rho_a = meanr-rng+2*rng*rndu(numpri,1);  
                                   /* Creates "numpri" AR1 coefficients in 
                                    each iteration.  The average AR1 
                                    coefficient is set to "meanr", with a 
                                    range of plus or minus "std".*/
rho_b = meanr;  

meanrho = meanrho|meanc(rho_a);


pp_a = zeros(numpri,1);            /* pp will store current observation of 
                                    each of the price series */
pp_b = zeros(numpri,1);            

p_a = 0;                        /* p will take the average over pp */
p_b = 0; 

i = 1;
do until i == sampsiz + 51;

err=sterr*rndn(numpri,1);   

@ intc=0;                   DGP1@
intc=sterr*rndn(numpri,1); @DGP2@

pp2_a=intc+rho_a.*pp_a+err;
                              /* Updating each AR1 series.  Errors are 
                                 assumed normal, mean zero, standard error
                                 equal to sterr. */
pp2_b=intc+rho_b.*pp_b+err;

x_a=meanc(pp2_a);                  /* Average of latest observations */
x_b=meanc(pp2_b);                  

p_a=p_a|x_a;                         /* stores the average price */
p_b=p_b|x_b;                         

pp_a = pp2_a;
pp_b = pp2_b;

i = i+1;
endo;

p_a=p_a[rows(p_a)-sampsiz:rows(p_a),1];   /* just take last "sampsiz" observations 
                                                                        to avoid any start-up bias */
p_b=p_b[rows(p_b)-sampsiz:rows(p_b),1];   

plag_a=p_a[1:rows(p_a)-1,1];
plag_b=p_b[1:rows(p_b)-1,1];

p_a=p_a[2:rows(p_a),1];
p_b=p_b[2:rows(p_b),1];

p_a=p_a-meanc(p_a);                 /* deviations from mean */
p_b=p_b-meanc(p_b);                 

plag_a=plag_a-meanc(plag_a);
plag_b=plag_b-meanc(plag_b);


slope_a=p_a/plag_a;                 /* estimates AR1 coefficient on the average
                                  price */

slope_b=p_b/plag_b;                 

coeffs_a=coeffs_a|slope_a;          /* store the slope estimate */
coeffs_b=coeffs_b|slope_b;          

j = j+1;
endo;

coeffs_a=coeffs_a[2:rows(coeffs_a),1];
coeffs_b=coeffs_b[2:rows(coeffs_b),1];

"Average AR1 coefficient estimated from average price: (Bias-uncorrected)";;meanc(coeffs_a);
@"Average AR1 coefficient estimated from average price: (no Aggregation Bias)";;meanc(coeffs_b);@

meanrho=meanrho[2:rows(meanrho),1];

"True average AR1 coefficient: ";; meanc(meanrho);

"Average AR1 coefficient estimated from average price: (Bias-corrected)";; meanc(coeffs_a)+(meanc(meanrho)-meanc(coeffs_b));