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.
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(235.1
+ ,1
+ ,280.7
+ ,1
+ ,264.6
+ ,2
+ ,240.7
+ ,0
+ ,201.4
+ ,1
+ ,240.8
+ ,0
+ ,241.1
+ ,-1
+ ,223.8
+ ,-3
+ ,206.1
+ ,-3
+ ,174.7
+ ,-3
+ ,203.3
+ ,-4
+ ,220.5
+ ,-8
+ ,299.5
+ ,-9
+ ,347.4
+ ,-13
+ ,338.3
+ ,-18
+ ,327.7
+ ,-11
+ ,351.6
+ ,-9
+ ,396.6
+ ,-10
+ ,438.8
+ ,-13
+ ,395.6
+ ,-11
+ ,363.5
+ ,-5
+ ,378.8
+ ,-15
+ ,357
+ ,-6
+ ,369
+ ,-6
+ ,464.8
+ ,-3
+ ,479.1
+ ,-1
+ ,431.3
+ ,-3
+ ,366.5
+ ,-4
+ ,326.3
+ ,-6
+ ,355.1
+ ,0
+ ,331.6
+ ,-4
+ ,261.3
+ ,-2
+ ,249
+ ,-2
+ ,205.5
+ ,-6
+ ,235.6
+ ,-7
+ ,240.9
+ ,-6
+ ,264.9
+ ,-6
+ ,253.8
+ ,-3
+ ,232.3
+ ,-2
+ ,193.8
+ ,-5
+ ,177
+ ,-11
+ ,213.2
+ ,-11
+ ,207.2
+ ,-11
+ ,180.6
+ ,-10
+ ,188.6
+ ,-14
+ ,175.4
+ ,-8
+ ,199
+ ,-9
+ ,179.6
+ ,-5
+ ,225.8
+ ,-1
+ ,234
+ ,-2
+ ,200.2
+ ,-5
+ ,183.6
+ ,-4
+ ,178.2
+ ,-6
+ ,203.2
+ ,-2
+ ,208.5
+ ,-2
+ ,191.8
+ ,-2
+ ,172.8
+ ,-2
+ ,148
+ ,2
+ ,159.4
+ ,1
+ ,154.5
+ ,-8
+ ,213.2
+ ,-1
+ ,196.4
+ ,1
+ ,182.8
+ ,-1
+ ,176.4
+ ,2
+ ,153.6
+ ,2
+ ,173.2
+ ,1
+ ,171
+ ,-1
+ ,151.2
+ ,-2
+ ,161.9
+ ,-2
+ ,157.2
+ ,-1
+ ,201.7
+ ,-8
+ ,236.4
+ ,-4
+ ,356.1
+ ,-6
+ ,398.3
+ ,-3
+ ,403.7
+ ,-3
+ ,384.6
+ ,-7
+ ,365.8
+ ,-9
+ ,368.1
+ ,-11
+ ,367.9
+ ,-13
+ ,347
+ ,-11
+ ,343.3
+ ,-9
+ ,292.9
+ ,-17
+ ,311.5
+ ,-22
+ ,300.9
+ ,-25
+ ,366.9
+ ,-20
+ ,356.9
+ ,-24
+ ,329.7
+ ,-24
+ ,316.2
+ ,-22
+ ,269
+ ,-19
+ ,289.3
+ ,-18
+ ,266.2
+ ,-17
+ ,253.6
+ ,-11
+ ,233.8
+ ,-11
+ ,228.4
+ ,-12
+ ,253.6
+ ,-10
+ ,260.1
+ ,-15
+ ,306.6
+ ,-15
+ ,309.2
+ ,-15
+ ,309.5
+ ,-13
+ ,271
+ ,-8
+ ,279.9
+ ,-13
+ ,317.9
+ ,-9
+ ,298.4
+ ,-7
+ ,246.7
+ ,-4
+ ,227.3
+ ,-4
+ ,209.1
+ ,-2)
+ ,dim=c(2
+ ,106)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:106))
> y <- array(NA,dim=c(2,106),dimnames=list(c('Y','X'),1:106))
> 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 = 'Do not include Seasonal 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> 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
1 235.1 1
2 280.7 1
3 264.6 2
4 240.7 0
5 201.4 1
6 240.8 0
7 241.1 -1
8 223.8 -3
9 206.1 -3
10 174.7 -3
11 203.3 -4
12 220.5 -8
13 299.5 -9
14 347.4 -13
15 338.3 -18
16 327.7 -11
17 351.6 -9
18 396.6 -10
19 438.8 -13
20 395.6 -11
21 363.5 -5
22 378.8 -15
23 357.0 -6
24 369.0 -6
25 464.8 -3
26 479.1 -1
27 431.3 -3
28 366.5 -4
29 326.3 -6
30 355.1 0
31 331.6 -4
32 261.3 -2
33 249.0 -2
34 205.5 -6
35 235.6 -7
36 240.9 -6
37 264.9 -6
38 253.8 -3
39 232.3 -2
40 193.8 -5
41 177.0 -11
42 213.2 -11
43 207.2 -11
44 180.6 -10
45 188.6 -14
46 175.4 -8
47 199.0 -9
48 179.6 -5
49 225.8 -1
50 234.0 -2
51 200.2 -5
52 183.6 -4
53 178.2 -6
54 203.2 -2
55 208.5 -2
56 191.8 -2
57 172.8 -2
58 148.0 2
59 159.4 1
60 154.5 -8
61 213.2 -1
62 196.4 1
63 182.8 -1
64 176.4 2
65 153.6 2
66 173.2 1
67 171.0 -1
68 151.2 -2
69 161.9 -2
70 157.2 -1
71 201.7 -8
72 236.4 -4
73 356.1 -6
74 398.3 -3
75 403.7 -3
76 384.6 -7
77 365.8 -9
78 368.1 -11
79 367.9 -13
80 347.0 -11
81 343.3 -9
82 292.9 -17
83 311.5 -22
84 300.9 -25
85 366.9 -20
86 356.9 -24
87 329.7 -24
88 316.2 -22
89 269.0 -19
90 289.3 -18
91 266.2 -17
92 253.6 -11
93 233.8 -11
94 228.4 -12
95 253.6 -10
96 260.1 -15
97 306.6 -15
98 309.2 -15
99 309.5 -13
100 271.0 -8
101 279.9 -13
102 317.9 -9
103 298.4 -7
104 246.7 -4
105 227.3 -4
106 209.1 -2
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
235.632 -4.506
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-117.18 -56.09 -17.91 42.31 238.96
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 235.632 10.877 21.662 < 2e-16 ***
X -4.506 1.131 -3.984 0.000126 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 75.13 on 104 degrees of freedom
Multiple R-squared: 0.1324, Adjusted R-squared: 0.1241
F-statistic: 15.87 on 1 and 104 DF, p-value: 0.0001260
> 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.0984807878 1.969616e-01 9.015192e-01
[2,] 0.0349128709 6.982574e-02 9.650871e-01
[3,] 0.0117254142 2.345083e-02 9.882746e-01
[4,] 0.0034607745 6.921549e-03 9.965392e-01
[5,] 0.0011290405 2.258081e-03 9.988710e-01
[6,] 0.0008937363 1.787473e-03 9.991063e-01
[7,] 0.0002755430 5.510860e-04 9.997245e-01
[8,] 0.0003094902 6.189805e-04 9.996905e-01
[9,] 0.0031246496 6.249299e-03 9.968754e-01
[10,] 0.0082719731 1.654395e-02 9.917280e-01
[11,] 0.0043499834 8.699967e-03 9.956500e-01
[12,] 0.0029274669 5.854934e-03 9.970725e-01
[13,] 0.0040227226 8.045445e-03 9.959773e-01
[14,] 0.0121652298 2.433046e-02 9.878348e-01
[15,] 0.0330522247 6.610445e-02 9.669478e-01
[16,] 0.0371661755 7.433235e-02 9.628338e-01
[17,] 0.0525483646 1.050967e-01 9.474516e-01
[18,] 0.0377934347 7.558687e-02 9.622066e-01
[19,] 0.0407210713 8.144214e-02 9.592789e-01
[20,] 0.0500252350 1.000505e-01 9.499748e-01
[21,] 0.3687303619 7.374607e-01 6.312696e-01
[22,] 0.8589343250 2.821313e-01 1.410657e-01
[23,] 0.9550190702 8.996186e-02 4.498093e-02
[24,] 0.9647446764 7.051065e-02 3.525532e-02
[25,] 0.9588423549 8.231529e-02 4.115765e-02
[26,] 0.9745413536 5.091729e-02 2.545865e-02
[27,] 0.9746743870 5.065123e-02 2.532561e-02
[28,] 0.9683528367 6.329433e-02 3.164716e-02
[29,] 0.9611621705 7.767566e-02 3.883783e-02
[30,] 0.9686283189 6.274336e-02 3.137168e-02
[31,] 0.9673250700 6.534986e-02 3.267493e-02
[32,] 0.9628144339 7.437113e-02 3.718557e-02
[33,] 0.9541085263 9.178295e-02 4.589147e-02
[34,] 0.9431953722 1.136093e-01 5.680463e-02
[35,] 0.9314683699 1.370633e-01 6.853163e-02
[36,] 0.9386021325 1.227957e-01 6.139787e-02
[37,] 0.9693900248 6.121995e-02 3.060998e-02
[38,] 0.9738839147 5.223217e-02 2.611609e-02
[39,] 0.9780072006 4.398560e-02 2.199280e-02
[40,] 0.9853089454 2.938211e-02 1.469105e-02
[41,] 0.9912827428 1.743451e-02 8.717257e-03
[42,] 0.9935792091 1.284158e-02 6.420791e-03
[43,] 0.9938187457 1.236251e-02 6.181254e-03
[44,] 0.9942056837 1.158863e-02 5.794316e-03
[45,] 0.9919441608 1.611168e-02 8.055839e-03
[46,] 0.9887860683 2.242786e-02 1.121393e-02
[47,] 0.9869530118 2.609398e-02 1.304699e-02
[48,] 0.9862071952 2.758561e-02 1.379280e-02
[49,] 0.9872283611 2.554328e-02 1.277164e-02
[50,] 0.9835211924 3.295762e-02 1.647881e-02
[51,] 0.9783632992 4.327340e-02 2.163670e-02
[52,] 0.9738500772 5.229985e-02 2.614992e-02
[53,] 0.9722112571 5.557749e-02 2.778874e-02
[54,] 0.9718467110 5.630658e-02 2.815329e-02
[55,] 0.9698003670 6.039927e-02 3.019963e-02
[56,] 0.9815229235 3.695415e-02 1.847708e-02
[57,] 0.9746962357 5.060753e-02 2.530376e-02
[58,] 0.9666486515 6.670270e-02 3.335135e-02
[59,] 0.9614119169 7.717617e-02 3.858808e-02
[60,] 0.9539210325 9.215793e-02 4.607897e-02
[61,] 0.9545597834 9.088043e-02 4.544022e-02
[62,] 0.9512668561 9.746629e-02 4.873314e-02
[63,] 0.9542058832 9.158823e-02 4.579412e-02
[64,] 0.9704315787 5.913684e-02 2.956842e-02
[65,] 0.9814462283 3.710754e-02 1.855377e-02
[66,] 0.9917867110 1.642658e-02 8.213289e-03
[67,] 0.9947131851 1.057363e-02 5.286815e-03
[68,] 0.9942976769 1.140465e-02 5.702323e-03
[69,] 0.9940403268 1.191935e-02 5.959673e-03
[70,] 0.9978104556 4.379089e-03 2.189544e-03
[71,] 0.9996267085 7.465830e-04 3.732915e-04
[72,] 0.9999061356 1.877288e-04 9.386440e-05
[73,] 0.9999604511 7.909774e-05 3.954887e-05
[74,] 0.9999859379 2.812425e-05 1.406213e-05
[75,] 0.9999952999 9.400113e-06 4.700056e-06
[76,] 0.9999979552 4.089545e-06 2.044773e-06
[77,] 0.9999996260 7.480670e-07 3.740335e-07
[78,] 0.9999989704 2.059256e-06 1.029628e-06
[79,] 0.9999972924 5.415195e-06 2.707598e-06
[80,] 0.9999956755 8.648957e-06 4.324479e-06
[81,] 0.9999973905 5.219038e-06 2.609519e-06
[82,] 0.9999962113 7.577363e-06 3.788682e-06
[83,] 0.9999903213 1.935736e-05 9.678680e-06
[84,] 0.9999747258 5.054835e-05 2.527417e-05
[85,] 0.9999496388 1.007225e-04 5.036124e-05
[86,] 0.9998627620 2.744760e-04 1.372380e-04
[87,] 0.9997435417 5.129165e-04 2.564583e-04
[88,] 0.9994136894 1.172621e-03 5.863106e-04
[89,] 0.9992542744 1.491451e-03 7.457256e-04
[90,] 0.9995705798 8.588404e-04 4.294202e-04
[91,] 0.9990561849 1.887630e-03 9.438151e-04
[92,] 0.9994604369 1.079126e-03 5.395631e-04
[93,] 0.9982868740 3.426252e-03 1.713126e-03
[94,] 0.9949708335 1.005833e-02 5.029167e-03
[95,] 0.9841007029 3.179859e-02 1.589930e-02
[96,] 0.9537380259 9.252395e-02 4.626197e-02
[97,] 0.9967958573 6.408285e-03 3.204143e-03
> postscript(file="/var/www/html/rcomp/tmp/1vsyh1291025095.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/html/rcomp/tmp/2o2x21291025095.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/html/rcomp/tmp/3o2x21291025095.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/html/rcomp/tmp/4o2x21291025095.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/html/rcomp/tmp/5ztfn1291025095.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 = 106
Frequency = 1
1 2 3 4 5 6
3.9740960 49.5740960 37.9799593 5.0682327 -29.7259040 5.1682327
7 8 9 10 11 12
0.9623694 -25.3493573 -43.0493573 -74.4493573 -50.3552206 -51.1786738
13 14 15 16 17 18
23.3154629 53.1920096 21.5626931 42.5037363 75.4154629 115.9095996
19 20 21 22 23 24
144.5920096 110.4037363 105.3389161 75.5802830 94.3330528 106.3330528
25 26 27 28 29 30
215.6506427 238.9623694 182.1506427 112.8447794 63.6330528 119.4682327
31 32 33 34 35 36
77.9447794 16.6565060 4.3565060 -57.1669472 -31.5728105 -21.7669472
37 38 39 40 41 42
2.2330528 4.6506427 -12.3434940 -64.3610839 -108.1962637 -71.9962637
43 44 45 46 47 48
-77.9962637 -100.0904004 -110.1138537 -96.2786738 -77.1845371 -78.5610839
49 50 51 52 53 54
-14.3376306 -10.6434940 -57.9610839 -70.0552206 -84.4669472 -41.4434940
55 56 57 58 59 60
-36.1434940 -52.8434940 -71.8434940 -78.6200407 -71.7259040 -117.1786738
61 62 63 64 65 66
-26.9376306 -34.7259040 -57.3376306 -50.2200407 -73.0200407 -57.9259040
67 68 69 70 71 72
-69.1376306 -93.4434940 -82.7434940 -82.9376306 -69.9786738 -17.2552206
73 74 75 76 77 78
93.4330528 149.1506427 154.5506427 117.4271895 89.6154629 82.9037363
79 80 81 82 83 84
73.6920096 61.8037363 67.1154629 -19.3314436 -23.2607602 -47.3783501
85 86 87 88 89 90
41.1509665 13.1275132 -14.0724868 -18.5607602 -52.2431702 -27.4373069
91 92 93 94 95 96
-46.0314436 -31.5962637 -51.3962637 -61.3021271 -27.0904004 -43.1197170
97 98 99 100 101 102
3.3802830 5.9802830 15.2920096 -0.6786738 -14.3079904 41.7154629
103 104 105 106
31.2271895 -6.9552206 -26.3552206 -35.5434940
> postscript(file="/var/www/html/rcomp/tmp/6ztfn1291025095.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 = 106
Frequency = 1
lag(myerror, k = 1) myerror
0 3.9740960 NA
1 49.5740960 3.9740960
2 37.9799593 49.5740960
3 5.0682327 37.9799593
4 -29.7259040 5.0682327
5 5.1682327 -29.7259040
6 0.9623694 5.1682327
7 -25.3493573 0.9623694
8 -43.0493573 -25.3493573
9 -74.4493573 -43.0493573
10 -50.3552206 -74.4493573
11 -51.1786738 -50.3552206
12 23.3154629 -51.1786738
13 53.1920096 23.3154629
14 21.5626931 53.1920096
15 42.5037363 21.5626931
16 75.4154629 42.5037363
17 115.9095996 75.4154629
18 144.5920096 115.9095996
19 110.4037363 144.5920096
20 105.3389161 110.4037363
21 75.5802830 105.3389161
22 94.3330528 75.5802830
23 106.3330528 94.3330528
24 215.6506427 106.3330528
25 238.9623694 215.6506427
26 182.1506427 238.9623694
27 112.8447794 182.1506427
28 63.6330528 112.8447794
29 119.4682327 63.6330528
30 77.9447794 119.4682327
31 16.6565060 77.9447794
32 4.3565060 16.6565060
33 -57.1669472 4.3565060
34 -31.5728105 -57.1669472
35 -21.7669472 -31.5728105
36 2.2330528 -21.7669472
37 4.6506427 2.2330528
38 -12.3434940 4.6506427
39 -64.3610839 -12.3434940
40 -108.1962637 -64.3610839
41 -71.9962637 -108.1962637
42 -77.9962637 -71.9962637
43 -100.0904004 -77.9962637
44 -110.1138537 -100.0904004
45 -96.2786738 -110.1138537
46 -77.1845371 -96.2786738
47 -78.5610839 -77.1845371
48 -14.3376306 -78.5610839
49 -10.6434940 -14.3376306
50 -57.9610839 -10.6434940
51 -70.0552206 -57.9610839
52 -84.4669472 -70.0552206
53 -41.4434940 -84.4669472
54 -36.1434940 -41.4434940
55 -52.8434940 -36.1434940
56 -71.8434940 -52.8434940
57 -78.6200407 -71.8434940
58 -71.7259040 -78.6200407
59 -117.1786738 -71.7259040
60 -26.9376306 -117.1786738
61 -34.7259040 -26.9376306
62 -57.3376306 -34.7259040
63 -50.2200407 -57.3376306
64 -73.0200407 -50.2200407
65 -57.9259040 -73.0200407
66 -69.1376306 -57.9259040
67 -93.4434940 -69.1376306
68 -82.7434940 -93.4434940
69 -82.9376306 -82.7434940
70 -69.9786738 -82.9376306
71 -17.2552206 -69.9786738
72 93.4330528 -17.2552206
73 149.1506427 93.4330528
74 154.5506427 149.1506427
75 117.4271895 154.5506427
76 89.6154629 117.4271895
77 82.9037363 89.6154629
78 73.6920096 82.9037363
79 61.8037363 73.6920096
80 67.1154629 61.8037363
81 -19.3314436 67.1154629
82 -23.2607602 -19.3314436
83 -47.3783501 -23.2607602
84 41.1509665 -47.3783501
85 13.1275132 41.1509665
86 -14.0724868 13.1275132
87 -18.5607602 -14.0724868
88 -52.2431702 -18.5607602
89 -27.4373069 -52.2431702
90 -46.0314436 -27.4373069
91 -31.5962637 -46.0314436
92 -51.3962637 -31.5962637
93 -61.3021271 -51.3962637
94 -27.0904004 -61.3021271
95 -43.1197170 -27.0904004
96 3.3802830 -43.1197170
97 5.9802830 3.3802830
98 15.2920096 5.9802830
99 -0.6786738 15.2920096
100 -14.3079904 -0.6786738
101 41.7154629 -14.3079904
102 31.2271895 41.7154629
103 -6.9552206 31.2271895
104 -26.3552206 -6.9552206
105 -35.5434940 -26.3552206
106 NA -35.5434940
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 49.5740960 3.9740960
[2,] 37.9799593 49.5740960
[3,] 5.0682327 37.9799593
[4,] -29.7259040 5.0682327
[5,] 5.1682327 -29.7259040
[6,] 0.9623694 5.1682327
[7,] -25.3493573 0.9623694
[8,] -43.0493573 -25.3493573
[9,] -74.4493573 -43.0493573
[10,] -50.3552206 -74.4493573
[11,] -51.1786738 -50.3552206
[12,] 23.3154629 -51.1786738
[13,] 53.1920096 23.3154629
[14,] 21.5626931 53.1920096
[15,] 42.5037363 21.5626931
[16,] 75.4154629 42.5037363
[17,] 115.9095996 75.4154629
[18,] 144.5920096 115.9095996
[19,] 110.4037363 144.5920096
[20,] 105.3389161 110.4037363
[21,] 75.5802830 105.3389161
[22,] 94.3330528 75.5802830
[23,] 106.3330528 94.3330528
[24,] 215.6506427 106.3330528
[25,] 238.9623694 215.6506427
[26,] 182.1506427 238.9623694
[27,] 112.8447794 182.1506427
[28,] 63.6330528 112.8447794
[29,] 119.4682327 63.6330528
[30,] 77.9447794 119.4682327
[31,] 16.6565060 77.9447794
[32,] 4.3565060 16.6565060
[33,] -57.1669472 4.3565060
[34,] -31.5728105 -57.1669472
[35,] -21.7669472 -31.5728105
[36,] 2.2330528 -21.7669472
[37,] 4.6506427 2.2330528
[38,] -12.3434940 4.6506427
[39,] -64.3610839 -12.3434940
[40,] -108.1962637 -64.3610839
[41,] -71.9962637 -108.1962637
[42,] -77.9962637 -71.9962637
[43,] -100.0904004 -77.9962637
[44,] -110.1138537 -100.0904004
[45,] -96.2786738 -110.1138537
[46,] -77.1845371 -96.2786738
[47,] -78.5610839 -77.1845371
[48,] -14.3376306 -78.5610839
[49,] -10.6434940 -14.3376306
[50,] -57.9610839 -10.6434940
[51,] -70.0552206 -57.9610839
[52,] -84.4669472 -70.0552206
[53,] -41.4434940 -84.4669472
[54,] -36.1434940 -41.4434940
[55,] -52.8434940 -36.1434940
[56,] -71.8434940 -52.8434940
[57,] -78.6200407 -71.8434940
[58,] -71.7259040 -78.6200407
[59,] -117.1786738 -71.7259040
[60,] -26.9376306 -117.1786738
[61,] -34.7259040 -26.9376306
[62,] -57.3376306 -34.7259040
[63,] -50.2200407 -57.3376306
[64,] -73.0200407 -50.2200407
[65,] -57.9259040 -73.0200407
[66,] -69.1376306 -57.9259040
[67,] -93.4434940 -69.1376306
[68,] -82.7434940 -93.4434940
[69,] -82.9376306 -82.7434940
[70,] -69.9786738 -82.9376306
[71,] -17.2552206 -69.9786738
[72,] 93.4330528 -17.2552206
[73,] 149.1506427 93.4330528
[74,] 154.5506427 149.1506427
[75,] 117.4271895 154.5506427
[76,] 89.6154629 117.4271895
[77,] 82.9037363 89.6154629
[78,] 73.6920096 82.9037363
[79,] 61.8037363 73.6920096
[80,] 67.1154629 61.8037363
[81,] -19.3314436 67.1154629
[82,] -23.2607602 -19.3314436
[83,] -47.3783501 -23.2607602
[84,] 41.1509665 -47.3783501
[85,] 13.1275132 41.1509665
[86,] -14.0724868 13.1275132
[87,] -18.5607602 -14.0724868
[88,] -52.2431702 -18.5607602
[89,] -27.4373069 -52.2431702
[90,] -46.0314436 -27.4373069
[91,] -31.5962637 -46.0314436
[92,] -51.3962637 -31.5962637
[93,] -61.3021271 -51.3962637
[94,] -27.0904004 -61.3021271
[95,] -43.1197170 -27.0904004
[96,] 3.3802830 -43.1197170
[97,] 5.9802830 3.3802830
[98,] 15.2920096 5.9802830
[99,] -0.6786738 15.2920096
[100,] -14.3079904 -0.6786738
[101,] 41.7154629 -14.3079904
[102,] 31.2271895 41.7154629
[103,] -6.9552206 31.2271895
[104,] -26.3552206 -6.9552206
[105,] -35.5434940 -26.3552206
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 49.5740960 3.9740960
2 37.9799593 49.5740960
3 5.0682327 37.9799593
4 -29.7259040 5.0682327
5 5.1682327 -29.7259040
6 0.9623694 5.1682327
7 -25.3493573 0.9623694
8 -43.0493573 -25.3493573
9 -74.4493573 -43.0493573
10 -50.3552206 -74.4493573
11 -51.1786738 -50.3552206
12 23.3154629 -51.1786738
13 53.1920096 23.3154629
14 21.5626931 53.1920096
15 42.5037363 21.5626931
16 75.4154629 42.5037363
17 115.9095996 75.4154629
18 144.5920096 115.9095996
19 110.4037363 144.5920096
20 105.3389161 110.4037363
21 75.5802830 105.3389161
22 94.3330528 75.5802830
23 106.3330528 94.3330528
24 215.6506427 106.3330528
25 238.9623694 215.6506427
26 182.1506427 238.9623694
27 112.8447794 182.1506427
28 63.6330528 112.8447794
29 119.4682327 63.6330528
30 77.9447794 119.4682327
31 16.6565060 77.9447794
32 4.3565060 16.6565060
33 -57.1669472 4.3565060
34 -31.5728105 -57.1669472
35 -21.7669472 -31.5728105
36 2.2330528 -21.7669472
37 4.6506427 2.2330528
38 -12.3434940 4.6506427
39 -64.3610839 -12.3434940
40 -108.1962637 -64.3610839
41 -71.9962637 -108.1962637
42 -77.9962637 -71.9962637
43 -100.0904004 -77.9962637
44 -110.1138537 -100.0904004
45 -96.2786738 -110.1138537
46 -77.1845371 -96.2786738
47 -78.5610839 -77.1845371
48 -14.3376306 -78.5610839
49 -10.6434940 -14.3376306
50 -57.9610839 -10.6434940
51 -70.0552206 -57.9610839
52 -84.4669472 -70.0552206
53 -41.4434940 -84.4669472
54 -36.1434940 -41.4434940
55 -52.8434940 -36.1434940
56 -71.8434940 -52.8434940
57 -78.6200407 -71.8434940
58 -71.7259040 -78.6200407
59 -117.1786738 -71.7259040
60 -26.9376306 -117.1786738
61 -34.7259040 -26.9376306
62 -57.3376306 -34.7259040
63 -50.2200407 -57.3376306
64 -73.0200407 -50.2200407
65 -57.9259040 -73.0200407
66 -69.1376306 -57.9259040
67 -93.4434940 -69.1376306
68 -82.7434940 -93.4434940
69 -82.9376306 -82.7434940
70 -69.9786738 -82.9376306
71 -17.2552206 -69.9786738
72 93.4330528 -17.2552206
73 149.1506427 93.4330528
74 154.5506427 149.1506427
75 117.4271895 154.5506427
76 89.6154629 117.4271895
77 82.9037363 89.6154629
78 73.6920096 82.9037363
79 61.8037363 73.6920096
80 67.1154629 61.8037363
81 -19.3314436 67.1154629
82 -23.2607602 -19.3314436
83 -47.3783501 -23.2607602
84 41.1509665 -47.3783501
85 13.1275132 41.1509665
86 -14.0724868 13.1275132
87 -18.5607602 -14.0724868
88 -52.2431702 -18.5607602
89 -27.4373069 -52.2431702
90 -46.0314436 -27.4373069
91 -31.5962637 -46.0314436
92 -51.3962637 -31.5962637
93 -61.3021271 -51.3962637
94 -27.0904004 -61.3021271
95 -43.1197170 -27.0904004
96 3.3802830 -43.1197170
97 5.9802830 3.3802830
98 15.2920096 5.9802830
99 -0.6786738 15.2920096
100 -14.3079904 -0.6786738
101 41.7154629 -14.3079904
102 31.2271895 41.7154629
103 -6.9552206 31.2271895
104 -26.3552206 -6.9552206
105 -35.5434940 -26.3552206
> 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/html/rcomp/tmp/7r2eq1291025095.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/html/rcomp/tmp/8r2eq1291025095.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/html/rcomp/tmp/92bvb1291025095.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/html/rcomp/tmp/102bvb1291025095.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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/html/rcomp/tmp/11nccz1291025095.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/html/rcomp/tmp/129cam1291025095.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/html/rcomp/tmp/13fwtq1291025096.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/html/rcomp/tmp/141frw1291025096.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/html/rcomp/tmp/15mxp21291025096.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/html/rcomp/tmp/167go81291025096.tab")
+ }
>
> try(system("convert tmp/1vsyh1291025095.ps tmp/1vsyh1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/2o2x21291025095.ps tmp/2o2x21291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/3o2x21291025095.ps tmp/3o2x21291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/4o2x21291025095.ps tmp/4o2x21291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ztfn1291025095.ps tmp/5ztfn1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ztfn1291025095.ps tmp/6ztfn1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/7r2eq1291025095.ps tmp/7r2eq1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/8r2eq1291025095.ps tmp/8r2eq1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/92bvb1291025095.ps tmp/92bvb1291025095.png",intern=TRUE))
character(0)
> try(system("convert tmp/102bvb1291025095.ps tmp/102bvb1291025095.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.065 1.674 6.884