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
Hi @DavoWW
Thank you so much for your assitance. I've tried this with my data, and it end well for me.
Thanks.
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.