rm(list = ls())
library(foreign)
library(car)
library(gplots)
install.packages('plm')
Error in install.packages : Updating loaded packages
library(plm)
data("Grunfeld", package = "plm")
Panel <- pdata.frame(Grunfeld, c("firm","year")) # set panel structure
summary(Panel)
      firm         year          inv              value            capital       
 1      :20   1935   : 10   Min.   :   0.93   Min.   :  58.12   Min.   :   0.80  
 2      :20   1936   : 10   1st Qu.:  33.56   1st Qu.: 199.97   1st Qu.:  79.17  
 3      :20   1937   : 10   Median :  57.48   Median : 517.95   Median : 205.60  
 4      :20   1938   : 10   Mean   : 145.96   Mean   :1081.68   Mean   : 276.02  
 5      :20   1939   : 10   3rd Qu.: 138.04   3rd Qu.:1679.85   3rd Qu.: 358.10  
 6      :20   1940   : 10   Max.   :1486.70   Max.   :6241.70   Max.   :2226.30  
 (Other):80   (Other):140                                                        
# A panel of 220 observations (20 years for 11 firms)
# invest  - Gross investment in 1947 dollars
# value   - Market value as of Dec. 31 in 1947 dollars
# capital - Stock of plant and equipment in 1947 dollars
#1. Descriptive statistics
coplot(value ~ year|firm, type="b", data=Panel)        # Points and lines

# error bars show 95% confidence interval
plotmeans(value ~ firm, main="Heterogeneity across firms", data=Panel)

plotmeans(value ~ year, main="Heterogeneity across years", data=Panel)

#2. pooled model
ols<-lm(value~inv +capital, data = Panel)
summary(ols)

Call:
lm(formula = value ~ inv + capital, data = Panel)

Residuals:
     Min       1Q   Median       3Q      Max 
-2010.54  -339.50  -184.08    76.66  2707.84 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 410.8156    64.1419   6.405 1.08e-09 ***
inv           5.7598     0.2909  19.803  < 2e-16 ***
capital      -0.6153     0.2095  -2.937  0.00371 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 666.5 on 197 degrees of freedom
Multiple R-squared:  0.7455,    Adjusted R-squared:  0.7429 
F-statistic: 288.5 on 2 and 197 DF,  p-value: < 2.2e-16
#3. firm fixed effects
firm_fe <-lm(value~inv +capital + factor(firm)-1, data = Panel)
summary(firm_fe)

Call:
lm(formula = value ~ inv + capital + factor(firm) - 1, data = Panel)

Residuals:
    Min      1Q  Median      3Q     Max 
-807.78  -88.56   -7.23   76.21 1367.48 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
inv               2.8562     0.3075   9.288  < 2e-16 ***
capital          -0.5079     0.1404  -3.618 0.000381 ***
factor(firm)1  2926.5580   138.5974  21.116  < 2e-16 ***
factor(firm)2   949.1876   113.0063   8.399 1.08e-14 ***
factor(firm)3  1852.3960    69.8935  26.503  < 2e-16 ***
factor(firm)4   508.8034    62.2072   8.179 4.18e-14 ***
factor(firm)5   302.1638    80.8145   3.739 0.000245 ***
factor(firm)6   314.5649    60.9123   5.164 6.12e-07 ***
factor(firm)7   173.7996    68.7628   2.528 0.012310 *  
factor(firm)8   591.8985    60.5953   9.768  < 2e-16 ***
factor(firm)9   365.3017    68.1931   5.357 2.46e-07 ***
factor(firm)10   65.1287    60.0932   1.084 0.279844    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 268.7 on 188 degrees of freedom
Multiple R-squared:  0.9765,    Adjusted R-squared:  0.975 
F-statistic: 651.1 on 12 and 188 DF,  p-value: < 2.2e-16
# F test for significance of individual fixed effects
pFtest(value~inv +capital, data = Panel,
       effect = "individual", model = "within")

    F test for individual effects

data:  value ~ inv + capital
F = 113.76, df1 = 9, df2 = 188, p-value < 2.2e-16
alternative hypothesis: significant effects
# year fixed effects
year_fe <-lm(value~inv +capital + factor(year)-1, data = Panel)
summary(year_fe)

Call:
lm(formula = value ~ inv + capital + factor(year) - 1, data = Panel)

Residuals:
    Min      1Q  Median      3Q     Max 
-2190.9  -343.4  -142.6   112.4  2431.3 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
inv                5.6215     0.3047  18.448  < 2e-16 ***
capital           -0.2984     0.2505  -1.191  0.23511    
factor(year)1935 317.1238   215.8700   1.469  0.14358    
factor(year)1936 530.5330   216.4912   2.451  0.01523 *  
factor(year)1937 694.4634   217.0610   3.199  0.00163 ** 
factor(year)1938 452.7692   216.8299   2.088  0.03821 *  
factor(year)1939 673.9940   217.2237   3.103  0.00223 ** 
factor(year)1940 544.6940   217.4715   2.505  0.01316 *  
factor(year)1941 355.2886   218.2272   1.628  0.10528    
factor(year)1942 252.1110   218.6843   1.153  0.25052    
factor(year)1943 396.7571   218.9853   1.812  0.07170 .  
factor(year)1944 411.2561   218.9368   1.878  0.06196 .  
factor(year)1945 507.5727   219.1908   2.316  0.02172 *  
factor(year)1946 366.6214   219.9844   1.667  0.09736 .  
factor(year)1947 192.0668   223.0447   0.861  0.39033    
factor(year)1948 135.7126   225.3689   0.602  0.54782    
factor(year)1949 249.4954   228.1197   1.094  0.27556    
factor(year)1950 248.5819   229.3393   1.084  0.27987    
factor(year)1951 213.4401   230.0569   0.928  0.35478    
factor(year)1952 139.9284   234.3183   0.597  0.55115    
factor(year)1953  97.1670   240.4837   0.404  0.68666    
factor(year)1954  91.3851   248.4663   0.368  0.71346    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 680.6 on 178 degrees of freedom
Multiple R-squared:  0.8573,    Adjusted R-squared:  0.8397 
F-statistic: 48.62 on 22 and 178 DF,  p-value: < 2.2e-16
# F test for significance of time fixed effects
pFtest(value~inv +capital, data = Panel,
       effect = "time", model = "within")

    F test for time effects

data:  value ~ inv + capital
F = 0.57595, df1 = 19, df2 = 178, p-value = 0.9197
alternative hypothesis: significant effects
# year and time fixed effects
firm_year_fe<-lm(value~inv +capital +factor(firm)+ factor(year)-1, data = Panel)
summary(firm_year_fe)

Call:
lm(formula = value ~ inv + capital + factor(firm) + factor(year) - 
    1, data = Panel)

Residuals:
    Min      1Q  Median      3Q     Max 
-760.84 -105.29    4.94  129.14 1042.21 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
inv                 2.5694     0.3002   8.560 6.65e-15 ***
capital            -0.5885     0.1605  -3.666 0.000329 ***
factor(firm)1    2841.3257   147.5440  19.257  < 2e-16 ***
factor(firm)2     778.8056   129.7848   6.001 1.16e-08 ***
factor(firm)3    1602.1249    95.0532  16.855  < 2e-16 ***
factor(firm)4     231.4143    93.7647   2.468 0.014581 *  
factor(firm)5      47.2631   103.1136   0.458 0.647283    
factor(firm)6      27.0015    93.0555   0.290 0.772045    
factor(firm)7     -99.0248    94.7843  -1.045 0.297636    
factor(firm)8     299.2420    93.1847   3.211 0.001582 ** 
factor(firm)9      89.4670    94.5261   0.946 0.345255    
factor(firm)10   -245.3673    94.5322  -2.596 0.010274 *  
factor(year)1936  304.9451   108.3165   2.815 0.005452 ** 
factor(year)1937  540.9025   108.6010   4.981 1.55e-06 ***
factor(year)1938  172.6836   108.6597   1.589 0.113881    
factor(year)1939  407.6552   108.8791   3.744 0.000248 ***
factor(year)1940  379.2306   108.5158   3.495 0.000606 ***
factor(year)1941  275.7897   108.8553   2.534 0.012201 *  
factor(year)1942  128.3897   109.0554   1.177 0.240736    
factor(year)1943  260.9976   109.2877   2.388 0.018035 *  
factor(year)1944  284.5931   109.2154   2.606 0.009984 ** 
factor(year)1945  392.8830   109.3179   3.594 0.000427 ***
factor(year)1946  369.6591   109.5954   3.373 0.000922 ***
factor(year)1947  173.2198   111.2191   1.557 0.121231    
factor(year)1948  149.5766   112.4949   1.330 0.185432    
factor(year)1949  228.7342   114.6859   1.994 0.047712 *  
factor(year)1950  269.7332   115.1108   2.343 0.020281 *  
factor(year)1951  389.0969   114.4665   3.399 0.000843 ***
factor(year)1952  407.3494   116.5492   3.495 0.000605 ***
factor(year)1953  544.8882   119.5505   4.558 9.88e-06 ***
factor(year)1954  556.8600   123.8315   4.497 1.28e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 241.7 on 169 degrees of freedom
Multiple R-squared:  0.9829,    Adjusted R-squared:  0.9798 
F-statistic: 313.7 on 31 and 169 DF,  p-value: < 2.2e-16
# F test for significance of time and individual fixed effects
pFtest(value~inv +capital, data = Panel,
       effect = "twoways", model = "within")

    F test for twoways effects

data:  value ~ inv + capital
F = 47.486, df1 = 28, df2 = 169, p-value < 2.2e-16
alternative hypothesis: significant effects
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KCnJtKGxpc3QgPSBscygpKQpsaWJyYXJ5KGZvcmVpZ24pCmxpYnJhcnkoY2FyKQpsaWJyYXJ5KGdwbG90cykKaW5zdGFsbC5wYWNrYWdlcygncGxtJykKbGlicmFyeShwbG0pCgpkYXRhKCJHcnVuZmVsZCIsIHBhY2thZ2UgPSAicGxtIikKUGFuZWwgPC0gcGRhdGEuZnJhbWUoR3J1bmZlbGQsIGMoImZpcm0iLCJ5ZWFyIikpICMgc2V0IHBhbmVsIHN0cnVjdHVyZQpzdW1tYXJ5KFBhbmVsKQoKIyBBIHBhbmVsIG9mIDIyMCBvYnNlcnZhdGlvbnMgKDIwIHllYXJzIGZvciAxMSBmaXJtcykKCiMgaW52ZXN0ICAtIEdyb3NzIGludmVzdG1lbnQgaW4gMTk0NyBkb2xsYXJzCiMgdmFsdWUgICAtIE1hcmtldCB2YWx1ZSBhcyBvZiBEZWMuIDMxIGluIDE5NDcgZG9sbGFycwojIGNhcGl0YWwgLSBTdG9jayBvZiBwbGFudCBhbmQgZXF1aXBtZW50IGluIDE5NDcgZG9sbGFycwoKIzEuIERlc2NyaXB0aXZlIHN0YXRpc3RpY3MKY29wbG90KHZhbHVlIH4geWVhcnxmaXJtLCB0eXBlPSJiIiwgZGF0YT1QYW5lbCkgICAgICAgICMgUG9pbnRzIGFuZCBsaW5lcwoKIyBlcnJvciBiYXJzIHNob3cgOTUlIGNvbmZpZGVuY2UgaW50ZXJ2YWwKCnBsb3RtZWFucyh2YWx1ZSB+IGZpcm0sIG1haW49IkhldGVyb2dlbmVpdHkgYWNyb3NzIGZpcm1zIiwgZGF0YT1QYW5lbCkKcGxvdG1lYW5zKHZhbHVlIH4geWVhciwgbWFpbj0iSGV0ZXJvZ2VuZWl0eSBhY3Jvc3MgeWVhcnMiLCBkYXRhPVBhbmVsKQoKIzIuIHBvb2xlZCBtb2RlbApvbHM8LWxtKHZhbHVlfmludiArY2FwaXRhbCwgZGF0YSA9IFBhbmVsKQpzdW1tYXJ5KG9scykKCiMzLiBmaXJtIGZpeGVkIGVmZmVjdHMKZmlybV9mZSA8LWxtKHZhbHVlfmludiArY2FwaXRhbCArIGZhY3RvcihmaXJtKS0xLCBkYXRhID0gUGFuZWwpCnN1bW1hcnkoZmlybV9mZSkKIyBGIHRlc3QgZm9yIHNpZ25pZmljYW5jZSBvZiBpbmRpdmlkdWFsIGZpeGVkIGVmZmVjdHMKcEZ0ZXN0KHZhbHVlfmludiArY2FwaXRhbCwgZGF0YSA9IFBhbmVsLAogICAgICAgZWZmZWN0ID0gImluZGl2aWR1YWwiLCBtb2RlbCA9ICJ3aXRoaW4iKQoKIyB5ZWFyIGZpeGVkIGVmZmVjdHMKeWVhcl9mZSA8LWxtKHZhbHVlfmludiArY2FwaXRhbCArIGZhY3Rvcih5ZWFyKS0xLCBkYXRhID0gUGFuZWwpCnN1bW1hcnkoeWVhcl9mZSkKIyBGIHRlc3QgZm9yIHNpZ25pZmljYW5jZSBvZiB0aW1lIGZpeGVkIGVmZmVjdHMKcEZ0ZXN0KHZhbHVlfmludiArY2FwaXRhbCwgZGF0YSA9IFBhbmVsLAogICAgICAgZWZmZWN0ID0gInRpbWUiLCBtb2RlbCA9ICJ3aXRoaW4iKQoKIyB5ZWFyIGFuZCB0aW1lIGZpeGVkIGVmZmVjdHMKZmlybV95ZWFyX2ZlPC1sbSh2YWx1ZX5pbnYgK2NhcGl0YWwgK2ZhY3RvcihmaXJtKSsgZmFjdG9yKHllYXIpLTEsIGRhdGEgPSBQYW5lbCkKc3VtbWFyeShmaXJtX3llYXJfZmUpCiMgRiB0ZXN0IGZvciBzaWduaWZpY2FuY2Ugb2YgdGltZSBhbmQgaW5kaXZpZHVhbCBmaXhlZCBlZmZlY3RzCnBGdGVzdCh2YWx1ZX5pbnYgK2NhcGl0YWwsIGRhdGEgPSBQYW5lbCwKICAgICAgIGVmZmVjdCA9ICJ0d293YXlzIiwgbW9kZWwgPSAid2l0aGluIikKCgoKYGBgCgo=