For a MonteCarlo study, we would like to have some simulated data that follows the process that the command xtregar, fe identifies:
\[Y = X\beta + c_i + \epsilon_{i,t}\] with \[c_i\] individual fixed-effects
and with \[\epsilon_{i,t+1} = \rho \epsilon_{i,t} + \eta_{i,t}\] (where \[\eta_{i,t}\] are i.i.d. and follow a normal law.)

This should be simple but when we generate the data with some given parameters, and we analyse the simulated data, we do not get these parameters back with xtregar, fe ! (more exactly, results are significantly different from the initial parameters) What could go wrong ?

Here are the built data and code. Thanks in advance for any hint!

Code:
clear all
        *to get easily a panel structure : nr year
        set obs 100
        gen nr=_n
        expand 19
        gen year = 2000
        sort nr
        bysort nr: replace year=year[_n-1]+1 if _n!=1
        sort nr year
        xtset nr year
        save initial0.dta, replace
    *We consider some arsbitrary parameters:
    use initial0.dta, clear
        scalar the_rho = 0.3
        scalar the_sigma_epsilon = 2
        scalar the_sigma_eta = the_sigma_epsilon * sqrt(1-the_rho*the_rho)
        scalar the_sigma_c_i = 0.9
        matrix the_m = (0,0,0)
        matrix the_sd = (the_sigma_eta,the_sigma_epsilon,the_sigma_c_i)
    * We try to simulate 100 times the corresponding process
    set seed 89
    gen rho_emp = 0
    gen sigma_e_emp = 0
    gen sigma_std_ci = 0    
    forvalues i = 1/100 {
        drawnorm eta epsilon0 c_i, means(the_m) sds(the_sd)
     bysort nr: gen epsilon= epsilon0 if _n==1
     bysort nr: replace epsilon=eta + the_rho * epsilon[_n-1] if _n>1
     bysort nr: replace c_i = c_i[1]
    gen y = 5 + c_i + epsilon
        xtregar y, fe
        replace rho_emp = e(rho_ar) if _n==`i'
        replace sigma_e_emp = e(sigma_e) if _n==`i'
        replace sigma_std_ci = e(sigma_u) if _n==`i'
        drop epsilon epsilon epsilon0 y eta c_i
    }
    //comparing the results
    keep if _n<=100
    keep rho_emp sigma_e_emp sigma_std_ci
    gen constant=1
    reg rho_emp constant
        test (_cons=0.3)
    reg sigma_e_emp constant
        test (_cons=2)
    reg sigma_std_ci constant
        test (_cons=0.9)
Additional note: the line
Code:
 scalar the_sigma_eta = the_sigma_epsilon * sqrt(1-the_rho*the_rho)
comes from the stationnarity of $\epsilon_{i,t}$.


The reference for xtregar, fe can be found here: https://www.stata.com/manuals13/xtxtregar.pdf