C:\Rmisc>C:\R-2.9.0\bin\Rterm --vanilla R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pdf("spiritsggr.pdf") > > spirits <- read.fwf("C:\\data\\spirits.dat", width=c(rep(8,5)), + col.names=c("year","consume","pop","inc","popinc")) > > spirits year consume pop inc popinc 1 1870 1.9565 1.7669 1.9176 1.0853 2 1871 1.9794 1.7766 1.9059 1.0728 3 1872 2.0120 1.7764 1.8798 1.0582 4 1873 2.0449 1.7942 1.8727 1.0438 5 1874 2.0561 1.8156 1.8984 1.0456 6 1875 2.0678 1.8083 1.9137 1.0583 7 1876 2.0561 1.8083 1.9176 1.0604 8 1877 2.0428 1.8067 1.9176 1.0614 9 1878 2.0290 1.8166 1.9420 1.0690 10 1879 1.9980 1.8041 1.9547 1.0835 11 1880 1.9884 1.8053 1.9379 1.0735 12 1881 1.9835 1.8242 1.9462 1.0669 13 1882 1.9773 1.8395 1.9504 1.0603 14 1883 1.9748 1.8464 1.9504 1.0563 15 1884 1.9629 1.8492 1.9723 1.0666 16 1885 1.9396 1.8668 2.0000 1.0714 17 1886 1.9309 1.8783 2.0097 1.0700 18 1887 1.9271 1.8914 2.0146 1.0651 19 1888 1.9239 1.9166 2.0146 1.0511 20 1889 1.9414 1.9363 2.0097 1.0379 21 1890 1.9685 1.9548 2.0097 1.0281 22 1891 1.9727 1.9453 2.0097 1.0331 23 1892 1.9736 1.9292 2.0048 1.0392 24 1893 1.9499 1.9209 2.0097 1.0462 25 1894 1.9432 1.9510 2.0296 1.0403 26 1895 1.9569 1.9776 2.0399 1.0315 27 1896 1.9647 1.9814 2.0399 1.0295 28 1897 1.9710 1.9819 2.0296 1.0241 29 1898 1.9719 1.9828 2.0146 1.0160 30 1899 1.9956 2.0076 2.0245 1.0084 31 1900 2.0000 2.0000 2.0000 1.0000 32 1901 1.9904 1.9939 2.0048 1.0055 33 1902 1.9752 1.9933 2.0048 1.0058 34 1903 1.9494 1.9797 2.0000 1.0103 35 1904 1.9332 1.9772 1.9952 1.0091 36 1905 1.9139 1.9924 1.9952 1.0014 37 1906 1.9091 2.0117 1.9905 0.9895 38 1907 1.9139 2.0204 1.9813 0.9806 39 1908 1.8886 2.0018 1.9905 0.9944 40 1909 1.7945 2.0038 1.9859 0.9911 41 1910 1.7644 2.0099 2.0518 1.0208 42 1911 1.7817 2.0174 2.0474 1.0149 43 1912 1.7784 2.0279 2.0341 1.0031 44 1913 1.7945 2.0359 2.0255 0.9949 45 1914 1.7888 2.0216 2.0341 1.0062 46 1915 1.8751 1.9896 1.9445 0.9773 47 1916 1.7853 1.9843 1.9939 1.0048 48 1917 1.6075 1.9764 2.2082 1.1173 49 1918 1.5185 1.9965 2.2700 1.1370 50 1919 1.6513 2.0652 2.2430 1.0861 51 1920 1.6247 2.0369 2.2567 1.1079 52 1921 1.5391 1.9723 2.2988 1.1655 53 1922 1.4922 1.9797 2.3723 1.1983 54 1923 1.4606 2.0136 2.4105 1.1971 55 1924 1.4551 2.0165 2.4081 1.1942 56 1925 1.4425 2.0213 2.4081 1.1914 57 1926 1.4023 2.0206 2.4367 1.2059 58 1927 1.3991 2.0563 2.4284 1.1810 59 1928 1.3798 2.0579 2.4310 1.1813 60 1929 1.3782 2.0649 2.4363 1.1799 61 1930 1.3366 2.0582 2.4552 1.1929 62 1931 1.3026 2.0517 2.4838 1.2106 63 1932 1.2592 2.0491 2.4958 1.2180 64 1933 1.2365 2.0766 2.5048 1.2062 65 1934 1.2549 2.0890 2.5017 1.1976 66 1935 1.2527 2.1059 2.4958 1.1851 67 1936 1.2763 2.1205 2.4838 1.1713 68 1937 1.2906 2.1205 2.4636 1.1618 69 1938 1.2721 2.1182 2.4580 1.1604 > allmat <- as.matrix(spirits) > allmat year consume pop inc popinc [1,] 1870 1.9565 1.7669 1.9176 1.0853 [2,] 1871 1.9794 1.7766 1.9059 1.0728 [3,] 1872 2.0120 1.7764 1.8798 1.0582 [4,] 1873 2.0449 1.7942 1.8727 1.0438 [5,] 1874 2.0561 1.8156 1.8984 1.0456 [6,] 1875 2.0678 1.8083 1.9137 1.0583 [7,] 1876 2.0561 1.8083 1.9176 1.0604 [8,] 1877 2.0428 1.8067 1.9176 1.0614 [9,] 1878 2.0290 1.8166 1.9420 1.0690 [10,] 1879 1.9980 1.8041 1.9547 1.0835 [11,] 1880 1.9884 1.8053 1.9379 1.0735 [12,] 1881 1.9835 1.8242 1.9462 1.0669 [13,] 1882 1.9773 1.8395 1.9504 1.0603 [14,] 1883 1.9748 1.8464 1.9504 1.0563 [15,] 1884 1.9629 1.8492 1.9723 1.0666 [16,] 1885 1.9396 1.8668 2.0000 1.0714 [17,] 1886 1.9309 1.8783 2.0097 1.0700 [18,] 1887 1.9271 1.8914 2.0146 1.0651 [19,] 1888 1.9239 1.9166 2.0146 1.0511 [20,] 1889 1.9414 1.9363 2.0097 1.0379 [21,] 1890 1.9685 1.9548 2.0097 1.0281 [22,] 1891 1.9727 1.9453 2.0097 1.0331 [23,] 1892 1.9736 1.9292 2.0048 1.0392 [24,] 1893 1.9499 1.9209 2.0097 1.0462 [25,] 1894 1.9432 1.9510 2.0296 1.0403 [26,] 1895 1.9569 1.9776 2.0399 1.0315 [27,] 1896 1.9647 1.9814 2.0399 1.0295 [28,] 1897 1.9710 1.9819 2.0296 1.0241 [29,] 1898 1.9719 1.9828 2.0146 1.0160 [30,] 1899 1.9956 2.0076 2.0245 1.0084 [31,] 1900 2.0000 2.0000 2.0000 1.0000 [32,] 1901 1.9904 1.9939 2.0048 1.0055 [33,] 1902 1.9752 1.9933 2.0048 1.0058 [34,] 1903 1.9494 1.9797 2.0000 1.0103 [35,] 1904 1.9332 1.9772 1.9952 1.0091 [36,] 1905 1.9139 1.9924 1.9952 1.0014 [37,] 1906 1.9091 2.0117 1.9905 0.9895 [38,] 1907 1.9139 2.0204 1.9813 0.9806 [39,] 1908 1.8886 2.0018 1.9905 0.9944 [40,] 1909 1.7945 2.0038 1.9859 0.9911 [41,] 1910 1.7644 2.0099 2.0518 1.0208 [42,] 1911 1.7817 2.0174 2.0474 1.0149 [43,] 1912 1.7784 2.0279 2.0341 1.0031 [44,] 1913 1.7945 2.0359 2.0255 0.9949 [45,] 1914 1.7888 2.0216 2.0341 1.0062 [46,] 1915 1.8751 1.9896 1.9445 0.9773 [47,] 1916 1.7853 1.9843 1.9939 1.0048 [48,] 1917 1.6075 1.9764 2.2082 1.1173 [49,] 1918 1.5185 1.9965 2.2700 1.1370 [50,] 1919 1.6513 2.0652 2.2430 1.0861 [51,] 1920 1.6247 2.0369 2.2567 1.1079 [52,] 1921 1.5391 1.9723 2.2988 1.1655 [53,] 1922 1.4922 1.9797 2.3723 1.1983 [54,] 1923 1.4606 2.0136 2.4105 1.1971 [55,] 1924 1.4551 2.0165 2.4081 1.1942 [56,] 1925 1.4425 2.0213 2.4081 1.1914 [57,] 1926 1.4023 2.0206 2.4367 1.2059 [58,] 1927 1.3991 2.0563 2.4284 1.1810 [59,] 1928 1.3798 2.0579 2.4310 1.1813 [60,] 1929 1.3782 2.0649 2.4363 1.1799 [61,] 1930 1.3366 2.0582 2.4552 1.1929 [62,] 1931 1.3026 2.0517 2.4838 1.2106 [63,] 1932 1.2592 2.0491 2.4958 1.2180 [64,] 1933 1.2365 2.0766 2.5048 1.2062 [65,] 1934 1.2549 2.0890 2.5017 1.1976 [66,] 1935 1.2527 2.1059 2.4958 1.1851 [67,] 1936 1.2763 2.1205 2.4838 1.1713 [68,] 1937 1.2906 2.1205 2.4636 1.1618 [69,] 1938 1.2721 2.1182 2.4580 1.1604 > nyear <- nrow(allmat) > > x0 <- matrix(rep(1,nyear),ncol=1) > > x <- cbind(x0,allmat[,3],allmat[,4]) > y <- allmat[,2] > > beta_ols <- solve(t(x) %*% x) %*% t(x) %*% y > yhat_ols <- x %*% beta_ols > e_ols <- y-yhat_ols > s2_ols <- (t(e_ols) %*% e_ols)/(nyear-ncol(x)) > varb_ols <- s2_ols[1,1]*solve(t(x) %*% x) > seb_ols <- sqrt(diag(varb_ols)) > > > print(cbind(beta_ols,seb_ols)) seb_ols [1,] 4.6117077 0.15261611 [2,] -0.1184552 0.10885040 [3,] -1.2317419 0.05024342 > > plot(e_ols,type='o') > > gam_0 <- 0 > gam_1 <- 0 > gam_2 <- 0 > gam_3 <- 0 > > for (t in 1:nyear) { + gam_0 <- gam_0 + (e_ols[t]*e_ols[t]) + } > gam_0 <- gam_0/nyear > > for (t in 1:(nyear-1)) { + gam_1 <- gam_1 + (e_ols[t]*e_ols[t+1]) + } > gam_1 <- gam_1/nyear > > for (t in 1:(nyear-2)) { + gam_2 <- gam_2 + (e_ols[t]*e_ols[t+2]) + } > gam_2 <- gam_2/nyear > > for (t in 1:(nyear-3)) { + gam_3 <- gam_3 + (e_ols[t]*e_ols[t+3]) + } > gam_3 <- gam_3/nyear > > > gam2_1 <- gam_0 > gam2_2 <- rbind(cbind(gam_0,gam_1),cbind(gam_1,gam_0)) > gam2_3 <- rbind(cbind(gam_0,gam_1,gam_2),cbind(gam_1,gam_0,gam_1),cbind(gam_2,gam_1,gam_0)) > > igam2_1 <- solve(gam2_1) > igam2_2 <- solve(gam2_2) > igam2_3 <- solve(gam2_3) > > gam1_1 <- gam_1 > gam1_2 <- rbind(gam_1,gam_2) > gam1_3 <- rbind(gam_1,gam_2,gam_3) > > rho_1 <- -igam2_1 %*% gam1_1 > rho_2 <- -igam2_2 %*% gam1_2 > rho_3 <- -igam2_3 %*% gam1_3 > > sigma2_1 <- gam_0 +(t(rho_1) %*% gam1_1) > sigma2_2 <- gam_0 +(t(rho_2) %*% gam1_2) > sigma2_3 <- gam_0 +(t(rho_3) %*% gam1_3) > > p_1 <- chol(igam2_1) > p_2 <- chol(igam2_2) > p_3 <- chol(igam2_3) > > t_1 <- matrix(rep(0,nyear*nyear),ncol=nyear) > t_2 <- matrix(rep(0,nyear*nyear),ncol=nyear) > t_3 <- matrix(rep(0,nyear*nyear),ncol=nyear) > > q_1 <- 1 > q_2 <- 2 > q_3 <- 3 > > for (row in 1:1) { + for (col in 1:1) { + t_1[row,col] <- sqrt(sigma2_1) * p_1[row,col] + }} > for (row in 2:nyear) { + t_1[row,row] <- 1 + prevcols <- 0 + for (col in (row-q_1):(row-1)) { + t_1[row,col] <- rho_1[(q_1-prevcols),1] + prevcols <- prevcols+1 + }} > > for (row in 1:2) { + for (col in 1:2) { + t_2[row,col] <- sqrt(sigma2_2) * p_2[row,col] + }} > for (row in 3:nyear) { + t_2[row,row] <- 1 + prevcols <- 0 + for (col in (row-q_2):(row-1)) { + t_2[row,col] <- rho_2[(q_2-prevcols),1] + prevcols <- prevcols+1 + }} > > for (row in 1:3) { + for (col in 1:3) { + t_3[row,col] <- sqrt(sigma2_3) * p_3[row,col] + }} > for (row in 4:nyear) { + t_3[row,row] <- 1 + prevcols <- 0 + for (col in (row-q_3):(row-1)) { + t_3[row,col] <- rho_3[(q_3-prevcols),1] + prevcols <- prevcols+1 + }} > > bgls_1 <- solve(t(x) %*% t(t_1) %*% t_1 %*% x) %*% t(x) %*% t(t_1) %*% t_1 %*% y > bgls_2 <- solve(t(x) %*% t(t_2) %*% t_2 %*% x) %*% t(x) %*% t(t_2) %*% t_2 %*% y > bgls_3 <- solve(t(x) %*% t(t_3) %*% t_3 %*% x) %*% t(x) %*% t(t_3) %*% t_3 %*% y > > s2gls_1 <- (t(y-x%*%bgls_1) %*% (t(t_1)%*%t_1) %*% (y-x%*%bgls_1))/(nyear-3-1) > s2gls_2 <- (t(y-x%*%bgls_2) %*% (t(t_2)%*%t_2) %*% (y-x%*%bgls_2))/(nyear-3-2) > s2gls_3 <- (t(y-x%*%bgls_3) %*% (t(t_3)%*%t_3) %*% (y-x%*%bgls_3))/(nyear-3-3) > > print (cbind(s2gls_1,s2gls_2,s2gls_3)) [,1] [,2] [,3] [1,] 0.0007104334 0.0007137001 0.0007244935 > > c_1 <- solve(t(x) %*% t(t_1) %*% t_1 %*% x) > c_2 <- solve(t(x) %*% t(t_2) %*% t_2 %*% x) > c_3 <- solve(t(x) %*% t(t_3) %*% t_3 %*% x) > > sbgls_1 <- sqrt(s2gls_1 * diag(c_1)) > sbgls_2 <- sqrt(s2gls_2 * diag(c_2)) > sbgls_3 <- sqrt(s2gls_3 * diag(c_3)) > > print (cbind(bgls_1,sbgls_1,bgls_2,sbgls_2,bgls_3,sbgls_3)) sbgls_1 sbgls_2 sbgls_3 [1,] 3.8169806 0.25947282 3.7866756 0.26146872 3.7939688 0.2605707 [2,] 0.2050046 0.14423814 0.2152456 0.14404774 0.2194104 0.1436347 [3,] -1.1609976 0.07473075 -1.1565219 0.07480849 -1.1634833 0.0745510 > > s2rho_1 <- sigma2_1/(nyear-3-1) > s2rho_2 <- sigma2_2/(nyear-3-2) > s2rho_3 <- sigma2_3/(nyear-3-3) > > serho_11 <- sqrt(s2rho_1 %*% igam2_1[1,1]) > serho_21 <- sqrt(s2rho_2 %*% igam2_2[1,1]) > serho_22 <- sqrt(s2rho_2 %*% igam2_2[2,2]) > serho_31 <- sqrt(s2rho_3 %*% igam2_3[1,1]) > serho_32 <- sqrt(s2rho_3 %*% igam2_3[2,2]) > serho_33 <- sqrt(s2rho_3 %*% igam2_3[3,3]) > > > serho_1 <- serho_11 > serho_2 <- rbind(serho_21,serho_22) > serho_3 <- rbind(serho_31,serho_32,serho_33) > > trho_11 <- rho_1[1,1]/serho_11 > trho_21 <- rho_2[1,1]/serho_21 > trho_22 <- rho_2[2,1]/serho_22 > trho_31 <- rho_3[1,1]/serho_31 > trho_32 <- rho_3[2,1]/serho_32 > trho_33 <- rho_3[3,1]/serho_33 > > trho_1 <- trho_11 > trho_2 <- rbind(trho_21,trho_22) > trho_3 <- rbind(trho_31,trho_32,trho_33) > > print (cbind(rho_1,serho_1,trho_1)) [,1] [,2] [,3] [1,] -0.8523311 0.06487048 -13.13897 > print (cbind(rho_2,serho_2,trho_2)) [,1] [,2] [,3] gam_0 -0.81912003 0.1249051 -6.5579405 gam_1 -0.03896502 0.1249051 -0.3119571 > print (cbind(rho_3,serho_3,trho_3)) [,1] [,2] [,3] gam_0 -0.82047208 0.1259123 -6.5162193 gam_1 -0.06738783 0.1626872 -0.4142171 gam_2 0.03469920 0.1259123 0.2755823 > > e_gls1 <- t_1 %*% (y-x%*%bgls_1) > e_gls2 <- t_2 %*% (y-x%*%bgls_2) > e_gls3 <- t_3 %*% (y-x%*%bgls_3) > > plot(e_gls1,type='o') > plot(e_gls2,type='o') > plot(e_gls3,type='o') > > dev.off() null device 1 > >