How to run Fixed Effect with autocorrelation AR(1)

Hello,
I am using a fixed effect model but I would like to use or consider an AR(1) process in my fixed effect model.
Please find here below my R code

fe_model <- plm(formula =sp ~ bcour+gvex+tradeop+grosssaving+gdp,data = DataN, model = "within",index = c("country","year")) summary(femodel)
I know there an PanelAR Package but i don't know how to use it for my fixed effect model.
Any assistance please?

Thanks in advance

Hi @yann,
Try this reproducible example:

library(plm)
library(panelAR)

# Using example data from the {rlm} package
data(Rehm)
head(Rehm)
#>   ccode year     NURR     gini  mean_ur selfemp cum_right tradeunion deficit
#> 1   USA 2001 52.16667 28.68078 4.769855 7.35761  417.2900       12.8  -0.632
#> 2   USA 2002 50.75000 25.68726 5.841096 7.24036  511.7300       12.6  -3.970
#> 3   USA 2003 54.16667 24.64738 6.053730 7.56592  606.4700       12.4  -4.966
#> 4   USA 2004 53.83333 24.21621 5.592452 7.57404  701.2100       12.0  -4.437
#> 5   CAN 2001 63.41667 23.03409 7.830982 9.88345  384.3836       30.4   0.658
#> 6   CAN 2002 62.91667 21.98265 8.257591 9.75807  384.3836       30.3  -0.094
#>   tradeopen gdp_growth
#> 1  24.63934   1.093330
#> 2  24.47346   1.827312
#> 3  24.65713   2.502588
#> 4  26.28587   3.584586
#> 5  80.57254   1.783798
#> 6  79.40430   2.924532
# help(Rehm)

# Simple fixed effects panel model
plm.out <- plm(NURR ~ gini + mean_ur + selfemp + cum_right + tradeunion +
                      deficit + tradeopen + gdp_growth,
               data = Rehm)
summary(plm.out)
#> Oneway (individual) effect Within Model
#> 
#> Call:
#> plm(formula = NURR ~ gini + mean_ur + selfemp + cum_right + tradeunion + 
#>     deficit + tradeopen + gdp_growth, data = Rehm)
#> 
#> Unbalanced Panel: n = 20, T = 1-4, N = 75
#> 
#> Residuals:
#>      Min.   1st Qu.    Median   3rd Qu.      Max. 
#> -4.348047 -0.752280  0.049648  0.643790  3.854271 
#> 
#> Coefficients:
#>              Estimate Std. Error t-value Pr(>|t|)  
#> gini        0.0264597  0.1456456  0.1817  0.85662  
#> mean_ur     0.1346773  0.4244673  0.3173  0.75243  
#> selfemp     0.4613859  0.5230737  0.8821  0.38223  
#> cum_right   0.0015867  0.0036951  0.4294  0.66958  
#> tradeunion -0.5117134  0.3871797 -1.3216  0.19268  
#> deficit    -0.4648856  0.2157942 -2.1543  0.03637 *
#> tradeopen  -0.1845863  0.1018647 -1.8121  0.07637 .
#> gdp_growth -0.0859647  0.2459948 -0.3495  0.72831  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Total Sum of Squares:    169.15
#> Residual Sum of Squares: 134.86
#> R-Squared:      0.2027
#> Adj. R-Squared: -0.25533
#> F-statistic: 1.49358 on 8 and 47 DF, p-value: 0.18518
r.squared(plm.out)
#> [1] 0.202696

# Overall autocorrelation
out <- panelAR(NURR ~ gini + mean_ur + selfemp + cum_right + tradeunion +
                      deficit + tradeopen + gdp_growth,
               data=Rehm,
               panelVar='ccode',
               timeVar='year',
               autoCorr='ar1',
               panelCorrMethod='pcse',
               rho.na.rm=TRUE,
               panel.weight='t-1',
               bound.rho=TRUE)
#> Panel-specific correlations bounded to [-1,1]

summary(out)
#> 
#> Panel Regression with AR(1) Prais-Winsten correction and panel-corrected standard errors
#> 
#> Unbalanced Panel Design:                                             
#>  Total obs.:       75 Avg obs. per panel 3.75
#>  Number of panels: 20 Max obs. per panel 4   
#>  Number of times:  4  Min obs. per panel 1   
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 98.849636   9.463231  10.446 1.27e-15 ***
#> gini        -1.517884   0.383107  -3.962 0.000185 ***
#> mean_ur      0.030108   0.226004   0.133 0.894426    
#> selfemp      0.286876   0.095881   2.992 0.003895 ** 
#> cum_right   -0.012181   0.002015  -6.044 7.76e-08 ***
#> tradeunion   0.029370   0.047242   0.622 0.536292    
#> deficit      0.497009   0.242515   2.049 0.044401 *  
#> tradeopen    0.060517   0.022469   2.693 0.008959 ** 
#> gdp_growth  -0.702535   0.504288  -1.393 0.168258    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-squared: 0.8504
#> Wald statistic: 163.8146, Pr(>Chisq(8)): 0
summary(out)$r2
#> [1] 0.8503908

# Panel-specific autocorrelation
out2 <- panelAR(NURR ~ gini + mean_ur + selfemp + cum_right + tradeunion +
                       deficit + tradeopen + gdp_growth,
                data=Rehm,
                panelVar='ccode',
                timeVar='year',
                autoCorr='psar1',
                panelCorrMethod='pcse',
                rho.na.rm=TRUE,
                panel.weight='t-1',
                bound.rho=TRUE)
#> Setting panel-specific correlation to 0 for at least one panel because unable to estimate autocorrelation.
#> Panel-specific correlations bounded to [-1,1]

summary(out2)
#> 
#> Panel Regression with AR(1) Prais-Winsten correction and panel-corrected standard errors
#> 
#> Unbalanced Panel Design:                                             
#>  Total obs.:       75 Avg obs. per panel 3.75
#>  Number of panels: 20 Max obs. per panel 4   
#>  Number of times:  4  Min obs. per panel 1   
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 110.789655   5.849267  18.941  < 2e-16 ***
#> gini         -1.789562   0.215858  -8.290 8.00e-12 ***
#> mean_ur      -0.622511   0.333244  -1.868  0.06620 .  
#> selfemp       0.388129   0.060097   6.458 1.47e-08 ***
#> cum_right    -0.016928   0.003713  -4.559 2.29e-05 ***
#> tradeunion   -0.033242   0.028079  -1.184  0.24071    
#> deficit       0.272098   0.258629   1.052  0.29660    
#> tradeopen     0.094478   0.032690   2.890  0.00521 ** 
#> gdp_growth   -0.915797   0.417649  -2.193  0.03186 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-squared: 0.9847
#> Wald statistic: 153.8551, Pr(>Chisq(8)): 0
summary(out2)$r2
#> [1] 0.9846603

Created on 2020-12-11 by the reprex package (v0.3.0)

HTH

Hi,
Thanks so much for your help, this is exactly what i was looking for.
May i ask you one little question :
-First how can get hte specific fixed effect for all countries (in the ccode column) after PanelAr regression? because when I run for example a plm within model,i can get them with the below code:

> plm.out <- plm(NURR ~ gini + mean_ur + selfemp + cum_right + tradeunion +
                      deficit + tradeopen + gdp_growth,
               data = Rehm)
>fixef(plm.out)

Thanks

Hi @yann,
As you have discovered the fixef() method is not available for a {PanelAR} output. You will need to use the predict() function as shown in the following code (continuing from the example posted previously) and then do a bit of massaging of the output:

fixef(plm.out)
predict(plm.out)

fixef(out)  # No method for PanelAR output
#help("predict.panelAR")
predict(out2)
predict(out2, se.fit=TRUE, conf.interval=TRUE)

# Output is ordered alphabetically on the panelVar factor (in this case cccode)
out2.df <- predict(out2, se.fit=TRUE, conf.interval=TRUE)$fit

# Sort output
library(tidyverse)
out2.df %>%
  rownames_to_column(., var="row_num") %>%
  arrange(., as.numeric(row_num)) -> temp.df

# Merge fitted with raw data and draw graphs.
Rehm %>%
  bind_cols(., temp.df) %>%
  ggplot(., aes(x=year, y=fit)) +
    geom_point() +
    geom_line(col="blue") +
    geom_point(aes(x=year, y=NURR), col="red", shape="square") +
    facet_wrap(~ ccode)

HTH

1 Like

Hi @DavoWW
Thank you so much for your assitance. I've tried this with my data, and it end well for me.
Thanks.

Hi @DavoWW,

Thanks for your kind assistance,

Best Regards,

Yann

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.