我想执行一个估计程序，该程序使用给定部门中所有公司的数据，为期5年。

我可以轻松地循环执行此操作，但是由于估算过程需要相当长的时间，因此我想对其进行并行化。有什么办法吗？

我的数据如下：

     sale_log   cogs_log    ppegt_log     m_naics4 naics_2 gvkey year
1   3.9070198  2.5146032  3.192821715 9.290151e-02      72  1001 1983
2   4.1028774  2.7375141  3.517861329 1.067687e-01      72  1001 1984
3   4.5909863  3.2106595  3.975112703 2.511660e-01      72  1001 1985
4   3.2560391  2.7867256 -0.763368555 1.351031e-02      44  1003 1982
5   3.2966287  2.8088799 -0.305698649 1.151525e-02      44  1003 1983
6   3.2636907  2.8330357  0.154036559 8.699394e-03      44  1003 1984
7   3.7916480  3.2346849  0.887916936 1.351803e-02      44  1003 1985
8   4.1778028  3.5364473  1.177985972 1.761273e-02      44  1003 1986
9   4.1819066  3.7297111  1.393016951 1.686331e-02      44  1003 1987
10  4.0174411  3.6050022  1.479584215 1.601205e-02      44  1003 1988
11  3.4466429  2.9633579  1.312863013 8.888067e-03      44  1003 1989
12  3.0667367  2.6128805  0.909779173 2.102674e-02      42  1004 1965
13  3.2362968  2.8140391  1.430690273 2.050934e-02      42  1004 1966
14  3.1981990  2.8822097  1.721614365 1.702929e-02      42  1004 1967
15  3.9265031  3.6159280  2.399823853 2.559074e-02      42  1004 1968
16  4.3343438  4.0116068  2.592692585 3.649313e-02      42  1004 1969
17  4.5869564  4.3059855  2.772196529 4.743631e-02      42  1004 1970
18  4.7015486  4.3995561  2.875267240 5.155589e-02      42  1004 1971
19  5.0564414  4.7539697  3.218686385 6.863808e-02      42  1004 1972
20  5.4323873  5.1711531  3.350849771 8.272720e-02      42  1004 1973
21  5.2979696  5.0033437  3.383504340 6.726429e-02      42  1004 1974
22  5.3958779  5.1475985  3.475121024 1.534230e-01      42  1004 1975
23  5.5442635  5.3195666  3.517557041 1.674937e-01      42  1004 1976
24  5.6260795  5.3909462  3.694842501 1.711362e-01      42  1004 1977
25  5.8039766  5.5455887  3.895724689 1.836405e-01      42  1004 1978
26  5.8198831  5.5665980  3.960153940 1.700499e-01      42  1004 1979
27  5.7474447  5.4697019  3.943733263 1.520660e-01      42  1004 1980

其中gvkey是公司ID，naic是行业代码。

我写的代码：

theta=matrix(,60,23)
count=1

temp <- dat %>% select(
  "sale_log","cogs_log","ppegt_log","m_naics4","naics_2","gvkey","year"
)

for (i in 1960:2019) { # 5-year rolling sector-year specific production functions

  sub <- temp[between(temp$year,i-5,i),] # subset 5 years

  jcount <- 1

  for (j in sort(unique(sub$naics_2))) { # loop over sectors
    temp2 <- sub[sub$naics_2==j,]
    mdl <- prodestOP(
              Y=temp2$sale_log,fX=temp2$cogs_log,sX=temp2$ppegt_log,pX=temp2$cogs_log,cX=temp2$m_naics4,idvar=temp2$gvkey,timevar=temp2$year
              ) 
    theta[count,jcount] <- mdl@Model$FSbetas[2]

    jcount <- jcount+1
  }

  count <- count+1

}