R version 2.11.1 (2010-05-31)
Copyright (C) 2010 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.
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.
> x <- array(list(1579,0,2146,0,2462,0,3695,0,4831,0,5134,0,6250,0,5760,0,6249,0,2917,0,1741,0,2359,0,1511,1,2059,0,2635,0,2867,0,4403,0,5720,0,4502,0,5749,0,5627,0,2846,0,1762,0,2429,0,1169,0,2154,1,2249,0,2687,0,4359,0,5382,0,4459,0,6398,0,4596,0,3024,0,1887,0,2070,0,1351,0,2218,0,2461,1,3028,0,4784,0,4975,0,4607,0,6249,0,4809,0,3157,0,1910,0,2228,0,1594,0,2467,0,2222,0,3607,1,4685,0,4962,0,5770,0,5480,0,5000,0,3228,0,1993,0,2288,0,1580,0,2111,0,2192,0,3601,0,4665,1,4876,0,5813,0,5589,0,5331,0,3075,0,2002,0,2306,0,1507,0,1992,0,2487,0,3490,0,4647,0,5594,1,5611,0,5788,0,6204,0,3013,0,1931,0,2549,0,1504,0,2090,0,2702,0,2939,0,4500,0,6208,0,6415,1,5657,0,5964,0,3163,0,1997,0,2422,0,1376,0,2202,0,2683,0,3303,0,5202,0,5231,0,4880,0,7998,1,4977,0,3531,0,2025,0,2205,0,1442,0,2238,0,2179,0,3218,0,5139,0,4990,0,4914,0,6084,0,5672,1,3548,0,1793,0,2086,0),dim=c(2,120),dimnames=list(c('Y','X'),1:120))
> y <- array(NA,dim=c(2,120),dimnames=list(c('Y','X'),1:120))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1579 0 1 0 0 0 0 0 0 0 0 0 0
2 2146 0 0 1 0 0 0 0 0 0 0 0 0
3 2462 0 0 0 1 0 0 0 0 0 0 0 0
4 3695 0 0 0 0 1 0 0 0 0 0 0 0
5 4831 0 0 0 0 0 1 0 0 0 0 0 0
6 5134 0 0 0 0 0 0 1 0 0 0 0 0
7 6250 0 0 0 0 0 0 0 1 0 0 0 0
8 5760 0 0 0 0 0 0 0 0 1 0 0 0
9 6249 0 0 0 0 0 0 0 0 0 1 0 0
10 2917 0 0 0 0 0 0 0 0 0 0 1 0
11 1741 0 0 0 0 0 0 0 0 0 0 0 1
12 2359 0 0 0 0 0 0 0 0 0 0 0 0
13 1511 1 1 0 0 0 0 0 0 0 0 0 0
14 2059 0 0 1 0 0 0 0 0 0 0 0 0
15 2635 0 0 0 1 0 0 0 0 0 0 0 0
16 2867 0 0 0 0 1 0 0 0 0 0 0 0
17 4403 0 0 0 0 0 1 0 0 0 0 0 0
18 5720 0 0 0 0 0 0 1 0 0 0 0 0
19 4502 0 0 0 0 0 0 0 1 0 0 0 0
20 5749 0 0 0 0 0 0 0 0 1 0 0 0
21 5627 0 0 0 0 0 0 0 0 0 1 0 0
22 2846 0 0 0 0 0 0 0 0 0 0 1 0
23 1762 0 0 0 0 0 0 0 0 0 0 0 1
24 2429 0 0 0 0 0 0 0 0 0 0 0 0
25 1169 0 1 0 0 0 0 0 0 0 0 0 0
26 2154 1 0 1 0 0 0 0 0 0 0 0 0
27 2249 0 0 0 1 0 0 0 0 0 0 0 0
28 2687 0 0 0 0 1 0 0 0 0 0 0 0
29 4359 0 0 0 0 0 1 0 0 0 0 0 0
30 5382 0 0 0 0 0 0 1 0 0 0 0 0
31 4459 0 0 0 0 0 0 0 1 0 0 0 0
32 6398 0 0 0 0 0 0 0 0 1 0 0 0
33 4596 0 0 0 0 0 0 0 0 0 1 0 0
34 3024 0 0 0 0 0 0 0 0 0 0 1 0
35 1887 0 0 0 0 0 0 0 0 0 0 0 1
36 2070 0 0 0 0 0 0 0 0 0 0 0 0
37 1351 0 1 0 0 0 0 0 0 0 0 0 0
38 2218 0 0 1 0 0 0 0 0 0 0 0 0
39 2461 1 0 0 1 0 0 0 0 0 0 0 0
40 3028 0 0 0 0 1 0 0 0 0 0 0 0
41 4784 0 0 0 0 0 1 0 0 0 0 0 0
42 4975 0 0 0 0 0 0 1 0 0 0 0 0
43 4607 0 0 0 0 0 0 0 1 0 0 0 0
44 6249 0 0 0 0 0 0 0 0 1 0 0 0
45 4809 0 0 0 0 0 0 0 0 0 1 0 0
46 3157 0 0 0 0 0 0 0 0 0 0 1 0
47 1910 0 0 0 0 0 0 0 0 0 0 0 1
48 2228 0 0 0 0 0 0 0 0 0 0 0 0
49 1594 0 1 0 0 0 0 0 0 0 0 0 0
50 2467 0 0 1 0 0 0 0 0 0 0 0 0
51 2222 0 0 0 1 0 0 0 0 0 0 0 0
52 3607 1 0 0 0 1 0 0 0 0 0 0 0
53 4685 0 0 0 0 0 1 0 0 0 0 0 0
54 4962 0 0 0 0 0 0 1 0 0 0 0 0
55 5770 0 0 0 0 0 0 0 1 0 0 0 0
56 5480 0 0 0 0 0 0 0 0 1 0 0 0
57 5000 0 0 0 0 0 0 0 0 0 1 0 0
58 3228 0 0 0 0 0 0 0 0 0 0 1 0
59 1993 0 0 0 0 0 0 0 0 0 0 0 1
60 2288 0 0 0 0 0 0 0 0 0 0 0 0
61 1580 0 1 0 0 0 0 0 0 0 0 0 0
62 2111 0 0 1 0 0 0 0 0 0 0 0 0
63 2192 0 0 0 1 0 0 0 0 0 0 0 0
64 3601 0 0 0 0 1 0 0 0 0 0 0 0
65 4665 1 0 0 0 0 1 0 0 0 0 0 0
66 4876 0 0 0 0 0 0 1 0 0 0 0 0
67 5813 0 0 0 0 0 0 0 1 0 0 0 0
68 5589 0 0 0 0 0 0 0 0 1 0 0 0
69 5331 0 0 0 0 0 0 0 0 0 1 0 0
70 3075 0 0 0 0 0 0 0 0 0 0 1 0
71 2002 0 0 0 0 0 0 0 0 0 0 0 1
72 2306 0 0 0 0 0 0 0 0 0 0 0 0
73 1507 0 1 0 0 0 0 0 0 0 0 0 0
74 1992 0 0 1 0 0 0 0 0 0 0 0 0
75 2487 0 0 0 1 0 0 0 0 0 0 0 0
76 3490 0 0 0 0 1 0 0 0 0 0 0 0
77 4647 0 0 0 0 0 1 0 0 0 0 0 0
78 5594 1 0 0 0 0 0 1 0 0 0 0 0
79 5611 0 0 0 0 0 0 0 1 0 0 0 0
80 5788 0 0 0 0 0 0 0 0 1 0 0 0
81 6204 0 0 0 0 0 0 0 0 0 1 0 0
82 3013 0 0 0 0 0 0 0 0 0 0 1 0
83 1931 0 0 0 0 0 0 0 0 0 0 0 1
84 2549 0 0 0 0 0 0 0 0 0 0 0 0
85 1504 0 1 0 0 0 0 0 0 0 0 0 0
86 2090 0 0 1 0 0 0 0 0 0 0 0 0
87 2702 0 0 0 1 0 0 0 0 0 0 0 0
88 2939 0 0 0 0 1 0 0 0 0 0 0 0
89 4500 0 0 0 0 0 1 0 0 0 0 0 0
90 6208 0 0 0 0 0 0 1 0 0 0 0 0
91 6415 1 0 0 0 0 0 0 1 0 0 0 0
92 5657 0 0 0 0 0 0 0 0 1 0 0 0
93 5964 0 0 0 0 0 0 0 0 0 1 0 0
94 3163 0 0 0 0 0 0 0 0 0 0 1 0
95 1997 0 0 0 0 0 0 0 0 0 0 0 1
96 2422 0 0 0 0 0 0 0 0 0 0 0 0
97 1376 0 1 0 0 0 0 0 0 0 0 0 0
98 2202 0 0 1 0 0 0 0 0 0 0 0 0
99 2683 0 0 0 1 0 0 0 0 0 0 0 0
100 3303 0 0 0 0 1 0 0 0 0 0 0 0
101 5202 0 0 0 0 0 1 0 0 0 0 0 0
102 5231 0 0 0 0 0 0 1 0 0 0 0 0
103 4880 0 0 0 0 0 0 0 1 0 0 0 0
104 7998 1 0 0 0 0 0 0 0 1 0 0 0
105 4977 0 0 0 0 0 0 0 0 0 1 0 0
106 3531 0 0 0 0 0 0 0 0 0 0 1 0
107 2025 0 0 0 0 0 0 0 0 0 0 0 1
108 2205 0 0 0 0 0 0 0 0 0 0 0 0
109 1442 0 1 0 0 0 0 0 0 0 0 0 0
110 2238 0 0 1 0 0 0 0 0 0 0 0 0
111 2179 0 0 0 1 0 0 0 0 0 0 0 0
112 3218 0 0 0 0 1 0 0 0 0 0 0 0
113 5139 0 0 0 0 0 1 0 0 0 0 0 0
114 4990 0 0 0 0 0 0 1 0 0 0 0 0
115 4914 0 0 0 0 0 0 0 1 0 0 0 0
116 6084 0 0 0 0 0 0 0 0 1 0 0 0
117 5672 1 0 0 0 0 0 0 0 0 1 0 0
118 3548 0 0 0 0 0 0 0 0 0 0 1 0
119 1793 0 0 0 0 0 0 0 0 0 0 0 1
120 2086 0 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
2294.20 482.52 -881.15 -174.75 84.75 901.05
M5 M6 M7 M8 M9 M10
2379.05 2964.75 2979.65 3732.75 3100.45 856.00
M11
-390.10
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-814.85 -234.64 -12.77 140.54 1488.53
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2294.20 122.10 18.789 < 2e-16 ***
X 482.52 135.67 3.557 0.000561 ***
M1 -881.15 173.21 -5.087 1.56e-06 ***
M2 -174.75 173.21 -1.009 0.315308
M3 84.75 173.21 0.489 0.625652
M4 901.05 173.21 5.202 9.56e-07 ***
M5 2379.05 173.21 13.735 < 2e-16 ***
M6 2964.75 173.21 17.116 < 2e-16 ***
M7 2979.65 173.21 17.202 < 2e-16 ***
M8 3732.75 173.21 21.550 < 2e-16 ***
M9 3100.45 173.21 17.900 < 2e-16 ***
M10 856.00 172.68 4.957 2.70e-06 ***
M11 -390.10 172.68 -2.259 0.025906 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 386.1 on 107 degrees of freedom
Multiple R-squared: 0.9491, Adjusted R-squared: 0.9434
F-statistic: 166.3 on 12 and 107 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.5205002 0.95899954 0.47949977
[2,] 0.4574568 0.91491363 0.54254319
[3,] 0.4912862 0.98257234 0.50871383
[4,] 0.9745914 0.05081714 0.02540857
[5,] 0.9560078 0.08798435 0.04399218
[6,] 0.9528221 0.09435570 0.04717785
[7,] 0.9277334 0.14453312 0.07226656
[8,] 0.8912904 0.21741917 0.10870959
[9,] 0.8451397 0.30972050 0.15486025
[10,] 0.8160664 0.36786726 0.18393363
[11,] 0.7783920 0.44321604 0.22160802
[12,] 0.7332359 0.53352825 0.26676412
[13,] 0.7614994 0.47700120 0.23850060
[14,] 0.7204121 0.55917579 0.27958790
[15,] 0.6564239 0.68715223 0.34357612
[16,] 0.8155361 0.36892780 0.18446390
[17,] 0.8360083 0.32798342 0.16399171
[18,] 0.9642962 0.07140766 0.03570383
[19,] 0.9509746 0.09805074 0.04902537
[20,] 0.9324962 0.13500759 0.06750379
[21,] 0.9172717 0.16545665 0.08272833
[22,] 0.8903709 0.21925826 0.10962913
[23,] 0.8583561 0.28328771 0.14164386
[24,] 0.8448067 0.31038659 0.15519329
[25,] 0.8074397 0.38512067 0.19256034
[26,] 0.7710882 0.45782360 0.22891180
[27,] 0.7547267 0.49054663 0.24527332
[28,] 0.8033405 0.39331891 0.19665945
[29,] 0.7752965 0.44940691 0.22470346
[30,] 0.8242466 0.35150680 0.17575340
[31,] 0.7895434 0.42091318 0.21045659
[32,] 0.7449892 0.51002161 0.25501080
[33,] 0.6952485 0.60950292 0.30475146
[34,] 0.6491426 0.70171483 0.35085741
[35,] 0.6259526 0.74809470 0.37404735
[36,] 0.5838123 0.83237532 0.41618766
[37,] 0.5804191 0.83916176 0.41958088
[38,] 0.5229820 0.95403595 0.47701797
[39,] 0.4945997 0.98919943 0.50540029
[40,] 0.5686049 0.86279016 0.43139508
[41,] 0.6107125 0.77857494 0.38928747
[42,] 0.6043569 0.79128615 0.39564308
[43,] 0.5558873 0.88822540 0.44411270
[44,] 0.5019624 0.99607524 0.49803762
[45,] 0.4435765 0.88715301 0.55642349
[46,] 0.3928106 0.78562112 0.60718944
[47,] 0.3382604 0.67652071 0.66173964
[48,] 0.3039969 0.60799372 0.69600314
[49,] 0.2981933 0.59638653 0.70180674
[50,] 0.3926412 0.78528231 0.60735884
[51,] 0.3818301 0.76366030 0.61816985
[52,] 0.4525795 0.90515894 0.54742053
[53,] 0.4682253 0.93645059 0.53177470
[54,] 0.4143671 0.82873421 0.58563289
[55,] 0.3673784 0.73475672 0.63262164
[56,] 0.3156100 0.63122007 0.68438996
[57,] 0.2641488 0.52829752 0.73585124
[58,] 0.2180566 0.43611318 0.78194341
[59,] 0.1829454 0.36589080 0.81705460
[60,] 0.1469018 0.29380355 0.85309823
[61,] 0.1294991 0.25899826 0.87050087
[62,] 0.1070924 0.21418477 0.89290762
[63,] 0.1573935 0.31478697 0.84260651
[64,] 0.1753337 0.35066734 0.82466633
[65,] 0.1572773 0.31455464 0.84272268
[66,] 0.4094996 0.81899926 0.59050037
[67,] 0.3855467 0.77109336 0.61445332
[68,] 0.3238183 0.64763654 0.67618173
[69,] 0.2939061 0.58781225 0.70609388
[70,] 0.2400061 0.48001216 0.75999392
[71,] 0.1926824 0.38536471 0.80731764
[72,] 0.1666116 0.33322315 0.83338843
[73,] 0.1455365 0.29107303 0.85446348
[74,] 0.1683804 0.33676090 0.83161955
[75,] 0.4762935 0.95258698 0.52370651
[76,] 0.4960897 0.99217935 0.50391032
[77,] 0.6648117 0.67037662 0.33518831
[78,] 0.9800077 0.03998450 0.01999225
[79,] 0.9811080 0.03778396 0.01889198
[80,] 0.9674848 0.06503034 0.03251517
[81,] 0.9589422 0.08211559 0.04105779
[82,] 0.9303632 0.13927364 0.06963682
[83,] 0.8858167 0.22836656 0.11418328
[84,] 0.9064917 0.18701655 0.09350827
[85,] 0.8463548 0.30729044 0.15364522
[86,] 0.7706211 0.45875783 0.22937891
[87,] 0.6850258 0.62994838 0.31497419
[88,] 0.5495833 0.90083340 0.45041670
[89,] 0.9755529 0.04889418 0.02444709
> postscript(file="/var/www/rcomp/tmp/19uc31290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/29uc31290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/3klt61290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/4klt61290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/5klt61290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 120
Frequency = 1
1 2 3 4 5 6
165.951852 26.551852 83.051852 499.751852 157.751852 -124.948148
7 8 9 10 11 12
976.151852 -266.948148 854.351852 -233.200000 -163.100000 64.800000
13 14 15 16 17 18
-384.566667 -60.448148 256.051852 -328.248148 -270.248148 461.051852
19 20 21 22 23 24
-771.848148 -277.948148 232.351852 -304.200000 -142.100000 134.800000
25 26 27 28 29 30
-244.048148 -447.966667 -129.948148 -508.248148 -314.248148 123.051852
31 32 33 34 35 36
-814.848148 371.051852 -798.648148 -126.200000 -17.100000 -224.200000
37 38 39 40 41 42
-62.048148 98.551852 -400.466667 -167.248148 110.751852 -283.948148
43 44 45 46 47 48
-666.848148 222.051852 -585.648148 6.800000 5.900000 -66.200000
49 50 51 52 53 54
180.951852 347.551852 -156.948148 -70.766667 11.751852 -296.948148
55 56 57 58 59 60
496.151852 -546.948148 -394.648148 77.800000 88.900000 -6.200000
61 62 63 64 65 66
166.951852 -8.448148 -186.948148 405.751852 -490.766667 -382.948148
67 68 69 70 71 72
539.151852 -437.948148 -63.648148 -75.200000 97.900000 11.800000
73 74 75 76 77 78
93.951852 -127.448148 108.051852 294.751852 -26.248148 -147.466667
79 80 81 82 83 84
337.151852 -238.948148 809.351852 -137.200000 26.900000 254.800000
85 86 87 88 89 90
90.951852 -29.448148 323.051852 -256.248148 -173.248148 949.051852
91 92 93 94 95 96
658.633333 -369.948148 569.351852 12.800000 92.900000 127.800000
97 98 99 100 101 102
-37.048148 82.551852 304.051852 107.751852 528.751852 -27.948148
103 104 105 106 107 108
-393.848148 1488.533333 -417.648148 380.800000 120.900000 -89.200000
109 110 111 112 113 114
28.951852 118.551852 -199.948148 22.751852 465.751852 -268.948148
115 116 117 118 119 120
-359.848148 57.051852 -205.166667 397.800000 -111.100000 -208.200000
> postscript(file="/var/www/rcomp/tmp/6cusr1290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 120
Frequency = 1
lag(myerror, k = 1) myerror
0 165.951852 NA
1 26.551852 165.951852
2 83.051852 26.551852
3 499.751852 83.051852
4 157.751852 499.751852
5 -124.948148 157.751852
6 976.151852 -124.948148
7 -266.948148 976.151852
8 854.351852 -266.948148
9 -233.200000 854.351852
10 -163.100000 -233.200000
11 64.800000 -163.100000
12 -384.566667 64.800000
13 -60.448148 -384.566667
14 256.051852 -60.448148
15 -328.248148 256.051852
16 -270.248148 -328.248148
17 461.051852 -270.248148
18 -771.848148 461.051852
19 -277.948148 -771.848148
20 232.351852 -277.948148
21 -304.200000 232.351852
22 -142.100000 -304.200000
23 134.800000 -142.100000
24 -244.048148 134.800000
25 -447.966667 -244.048148
26 -129.948148 -447.966667
27 -508.248148 -129.948148
28 -314.248148 -508.248148
29 123.051852 -314.248148
30 -814.848148 123.051852
31 371.051852 -814.848148
32 -798.648148 371.051852
33 -126.200000 -798.648148
34 -17.100000 -126.200000
35 -224.200000 -17.100000
36 -62.048148 -224.200000
37 98.551852 -62.048148
38 -400.466667 98.551852
39 -167.248148 -400.466667
40 110.751852 -167.248148
41 -283.948148 110.751852
42 -666.848148 -283.948148
43 222.051852 -666.848148
44 -585.648148 222.051852
45 6.800000 -585.648148
46 5.900000 6.800000
47 -66.200000 5.900000
48 180.951852 -66.200000
49 347.551852 180.951852
50 -156.948148 347.551852
51 -70.766667 -156.948148
52 11.751852 -70.766667
53 -296.948148 11.751852
54 496.151852 -296.948148
55 -546.948148 496.151852
56 -394.648148 -546.948148
57 77.800000 -394.648148
58 88.900000 77.800000
59 -6.200000 88.900000
60 166.951852 -6.200000
61 -8.448148 166.951852
62 -186.948148 -8.448148
63 405.751852 -186.948148
64 -490.766667 405.751852
65 -382.948148 -490.766667
66 539.151852 -382.948148
67 -437.948148 539.151852
68 -63.648148 -437.948148
69 -75.200000 -63.648148
70 97.900000 -75.200000
71 11.800000 97.900000
72 93.951852 11.800000
73 -127.448148 93.951852
74 108.051852 -127.448148
75 294.751852 108.051852
76 -26.248148 294.751852
77 -147.466667 -26.248148
78 337.151852 -147.466667
79 -238.948148 337.151852
80 809.351852 -238.948148
81 -137.200000 809.351852
82 26.900000 -137.200000
83 254.800000 26.900000
84 90.951852 254.800000
85 -29.448148 90.951852
86 323.051852 -29.448148
87 -256.248148 323.051852
88 -173.248148 -256.248148
89 949.051852 -173.248148
90 658.633333 949.051852
91 -369.948148 658.633333
92 569.351852 -369.948148
93 12.800000 569.351852
94 92.900000 12.800000
95 127.800000 92.900000
96 -37.048148 127.800000
97 82.551852 -37.048148
98 304.051852 82.551852
99 107.751852 304.051852
100 528.751852 107.751852
101 -27.948148 528.751852
102 -393.848148 -27.948148
103 1488.533333 -393.848148
104 -417.648148 1488.533333
105 380.800000 -417.648148
106 120.900000 380.800000
107 -89.200000 120.900000
108 28.951852 -89.200000
109 118.551852 28.951852
110 -199.948148 118.551852
111 22.751852 -199.948148
112 465.751852 22.751852
113 -268.948148 465.751852
114 -359.848148 -268.948148
115 57.051852 -359.848148
116 -205.166667 57.051852
117 397.800000 -205.166667
118 -111.100000 397.800000
119 -208.200000 -111.100000
120 NA -208.200000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 26.551852 165.951852
[2,] 83.051852 26.551852
[3,] 499.751852 83.051852
[4,] 157.751852 499.751852
[5,] -124.948148 157.751852
[6,] 976.151852 -124.948148
[7,] -266.948148 976.151852
[8,] 854.351852 -266.948148
[9,] -233.200000 854.351852
[10,] -163.100000 -233.200000
[11,] 64.800000 -163.100000
[12,] -384.566667 64.800000
[13,] -60.448148 -384.566667
[14,] 256.051852 -60.448148
[15,] -328.248148 256.051852
[16,] -270.248148 -328.248148
[17,] 461.051852 -270.248148
[18,] -771.848148 461.051852
[19,] -277.948148 -771.848148
[20,] 232.351852 -277.948148
[21,] -304.200000 232.351852
[22,] -142.100000 -304.200000
[23,] 134.800000 -142.100000
[24,] -244.048148 134.800000
[25,] -447.966667 -244.048148
[26,] -129.948148 -447.966667
[27,] -508.248148 -129.948148
[28,] -314.248148 -508.248148
[29,] 123.051852 -314.248148
[30,] -814.848148 123.051852
[31,] 371.051852 -814.848148
[32,] -798.648148 371.051852
[33,] -126.200000 -798.648148
[34,] -17.100000 -126.200000
[35,] -224.200000 -17.100000
[36,] -62.048148 -224.200000
[37,] 98.551852 -62.048148
[38,] -400.466667 98.551852
[39,] -167.248148 -400.466667
[40,] 110.751852 -167.248148
[41,] -283.948148 110.751852
[42,] -666.848148 -283.948148
[43,] 222.051852 -666.848148
[44,] -585.648148 222.051852
[45,] 6.800000 -585.648148
[46,] 5.900000 6.800000
[47,] -66.200000 5.900000
[48,] 180.951852 -66.200000
[49,] 347.551852 180.951852
[50,] -156.948148 347.551852
[51,] -70.766667 -156.948148
[52,] 11.751852 -70.766667
[53,] -296.948148 11.751852
[54,] 496.151852 -296.948148
[55,] -546.948148 496.151852
[56,] -394.648148 -546.948148
[57,] 77.800000 -394.648148
[58,] 88.900000 77.800000
[59,] -6.200000 88.900000
[60,] 166.951852 -6.200000
[61,] -8.448148 166.951852
[62,] -186.948148 -8.448148
[63,] 405.751852 -186.948148
[64,] -490.766667 405.751852
[65,] -382.948148 -490.766667
[66,] 539.151852 -382.948148
[67,] -437.948148 539.151852
[68,] -63.648148 -437.948148
[69,] -75.200000 -63.648148
[70,] 97.900000 -75.200000
[71,] 11.800000 97.900000
[72,] 93.951852 11.800000
[73,] -127.448148 93.951852
[74,] 108.051852 -127.448148
[75,] 294.751852 108.051852
[76,] -26.248148 294.751852
[77,] -147.466667 -26.248148
[78,] 337.151852 -147.466667
[79,] -238.948148 337.151852
[80,] 809.351852 -238.948148
[81,] -137.200000 809.351852
[82,] 26.900000 -137.200000
[83,] 254.800000 26.900000
[84,] 90.951852 254.800000
[85,] -29.448148 90.951852
[86,] 323.051852 -29.448148
[87,] -256.248148 323.051852
[88,] -173.248148 -256.248148
[89,] 949.051852 -173.248148
[90,] 658.633333 949.051852
[91,] -369.948148 658.633333
[92,] 569.351852 -369.948148
[93,] 12.800000 569.351852
[94,] 92.900000 12.800000
[95,] 127.800000 92.900000
[96,] -37.048148 127.800000
[97,] 82.551852 -37.048148
[98,] 304.051852 82.551852
[99,] 107.751852 304.051852
[100,] 528.751852 107.751852
[101,] -27.948148 528.751852
[102,] -393.848148 -27.948148
[103,] 1488.533333 -393.848148
[104,] -417.648148 1488.533333
[105,] 380.800000 -417.648148
[106,] 120.900000 380.800000
[107,] -89.200000 120.900000
[108,] 28.951852 -89.200000
[109,] 118.551852 28.951852
[110,] -199.948148 118.551852
[111,] 22.751852 -199.948148
[112,] 465.751852 22.751852
[113,] -268.948148 465.751852
[114,] -359.848148 -268.948148
[115,] 57.051852 -359.848148
[116,] -205.166667 57.051852
[117,] 397.800000 -205.166667
[118,] -111.100000 397.800000
[119,] -208.200000 -111.100000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 26.551852 165.951852
2 83.051852 26.551852
3 499.751852 83.051852
4 157.751852 499.751852
5 -124.948148 157.751852
6 976.151852 -124.948148
7 -266.948148 976.151852
8 854.351852 -266.948148
9 -233.200000 854.351852
10 -163.100000 -233.200000
11 64.800000 -163.100000
12 -384.566667 64.800000
13 -60.448148 -384.566667
14 256.051852 -60.448148
15 -328.248148 256.051852
16 -270.248148 -328.248148
17 461.051852 -270.248148
18 -771.848148 461.051852
19 -277.948148 -771.848148
20 232.351852 -277.948148
21 -304.200000 232.351852
22 -142.100000 -304.200000
23 134.800000 -142.100000
24 -244.048148 134.800000
25 -447.966667 -244.048148
26 -129.948148 -447.966667
27 -508.248148 -129.948148
28 -314.248148 -508.248148
29 123.051852 -314.248148
30 -814.848148 123.051852
31 371.051852 -814.848148
32 -798.648148 371.051852
33 -126.200000 -798.648148
34 -17.100000 -126.200000
35 -224.200000 -17.100000
36 -62.048148 -224.200000
37 98.551852 -62.048148
38 -400.466667 98.551852
39 -167.248148 -400.466667
40 110.751852 -167.248148
41 -283.948148 110.751852
42 -666.848148 -283.948148
43 222.051852 -666.848148
44 -585.648148 222.051852
45 6.800000 -585.648148
46 5.900000 6.800000
47 -66.200000 5.900000
48 180.951852 -66.200000
49 347.551852 180.951852
50 -156.948148 347.551852
51 -70.766667 -156.948148
52 11.751852 -70.766667
53 -296.948148 11.751852
54 496.151852 -296.948148
55 -546.948148 496.151852
56 -394.648148 -546.948148
57 77.800000 -394.648148
58 88.900000 77.800000
59 -6.200000 88.900000
60 166.951852 -6.200000
61 -8.448148 166.951852
62 -186.948148 -8.448148
63 405.751852 -186.948148
64 -490.766667 405.751852
65 -382.948148 -490.766667
66 539.151852 -382.948148
67 -437.948148 539.151852
68 -63.648148 -437.948148
69 -75.200000 -63.648148
70 97.900000 -75.200000
71 11.800000 97.900000
72 93.951852 11.800000
73 -127.448148 93.951852
74 108.051852 -127.448148
75 294.751852 108.051852
76 -26.248148 294.751852
77 -147.466667 -26.248148
78 337.151852 -147.466667
79 -238.948148 337.151852
80 809.351852 -238.948148
81 -137.200000 809.351852
82 26.900000 -137.200000
83 254.800000 26.900000
84 90.951852 254.800000
85 -29.448148 90.951852
86 323.051852 -29.448148
87 -256.248148 323.051852
88 -173.248148 -256.248148
89 949.051852 -173.248148
90 658.633333 949.051852
91 -369.948148 658.633333
92 569.351852 -369.948148
93 12.800000 569.351852
94 92.900000 12.800000
95 127.800000 92.900000
96 -37.048148 127.800000
97 82.551852 -37.048148
98 304.051852 82.551852
99 107.751852 304.051852
100 528.751852 107.751852
101 -27.948148 528.751852
102 -393.848148 -27.948148
103 1488.533333 -393.848148
104 -417.648148 1488.533333
105 380.800000 -417.648148
106 120.900000 380.800000
107 -89.200000 120.900000
108 28.951852 -89.200000
109 118.551852 28.951852
110 -199.948148 118.551852
111 22.751852 -199.948148
112 465.751852 22.751852
113 -268.948148 465.751852
114 -359.848148 -268.948148
115 57.051852 -359.848148
116 -205.166667 57.051852
117 397.800000 -205.166667
118 -111.100000 397.800000
119 -208.200000 -111.100000
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/754sc1290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/854sc1290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/954sc1290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/10fd8f1290866485.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/11jdpl1290866485.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/12neor1290866485.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/131ol01290866485.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/14462n1290866485.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/157pib1290866485.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/16t7hh1290866485.tab")
+ }
>
> try(system("convert tmp/19uc31290866485.ps tmp/19uc31290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/29uc31290866485.ps tmp/29uc31290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/3klt61290866485.ps tmp/3klt61290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/4klt61290866485.ps tmp/4klt61290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/5klt61290866485.ps tmp/5klt61290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/6cusr1290866485.ps tmp/6cusr1290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/754sc1290866485.ps tmp/754sc1290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/854sc1290866485.ps tmp/854sc1290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/954sc1290866485.ps tmp/954sc1290866485.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fd8f1290866485.ps tmp/10fd8f1290866485.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.720 0.860 5.595