R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(1593,0,1477.9,0,1733.7,0,1569.7,0,1843.7,0,1950.3,0,1657.5,0,1772.1,0,1568.3,0,1809.8,0,1646.7,0,1808.5,0,1763.9,0,1625.5,0,1538.8,0,1342.4,0,1645.1,0,1619.9,0,1338.1,0,1505.5,0,1529.1,0,1511.9,0,1656.7,0,1694.4,0,1662.3,0,1588.7,0,1483.3,0,1585.6,0,1658.9,0,1584.4,0,1470.6,0,1618.7,0,1407.6,0,1473.9,0,1515.3,0,1485.4,0,1496.1,0,1493.5,0,1298.4,0,1375.3,0,1507.9,0,1455.3,0,1363.3,0,1392.8,0,1348.8,0,1880.3,0,1669.2,0,1543.6,0,1701.2,0,1516.5,0,1466.8,0,1484.1,0,1577.2,0,1684.5,0,1414.7,0,1674.5,0,1598.7,0,1739.1,0,1674.6,0,1671.8,0,1802,0,1526.8,0,1580.9,0,1634.8,0,1610.3,0,1712,0,1678.8,0,1708.1,0,1680.6,0,2056,1,1624,1,2021.4,1,1861.1,1,1750.8,1,1767.5,1,1710.3,1,2151.5,1,2047.9,1,1915.4,1,1984.7,1,1896.5,1,2170.8,1,2139.9,1,2330.5,1,2121.8,1,2226.8,1,1857.9,1,2155.9,1,2341.7,1,2290.2,1,2006.5,1,2111.9,1,1731.3,1,1762.2,1,1863.2,1,1943.5,1,1975.2,1),dim=c(2,97),dimnames=list(c('M','D'),1:97))
> y <- array(NA,dim=c(2,97),dimnames=list(c('M','D'),1:97))
> 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
M D
1 1593.0 0
2 1477.9 0
3 1733.7 0
4 1569.7 0
5 1843.7 0
6 1950.3 0
7 1657.5 0
8 1772.1 0
9 1568.3 0
10 1809.8 0
11 1646.7 0
12 1808.5 0
13 1763.9 0
14 1625.5 0
15 1538.8 0
16 1342.4 0
17 1645.1 0
18 1619.9 0
19 1338.1 0
20 1505.5 0
21 1529.1 0
22 1511.9 0
23 1656.7 0
24 1694.4 0
25 1662.3 0
26 1588.7 0
27 1483.3 0
28 1585.6 0
29 1658.9 0
30 1584.4 0
31 1470.6 0
32 1618.7 0
33 1407.6 0
34 1473.9 0
35 1515.3 0
36 1485.4 0
37 1496.1 0
38 1493.5 0
39 1298.4 0
40 1375.3 0
41 1507.9 0
42 1455.3 0
43 1363.3 0
44 1392.8 0
45 1348.8 0
46 1880.3 0
47 1669.2 0
48 1543.6 0
49 1701.2 0
50 1516.5 0
51 1466.8 0
52 1484.1 0
53 1577.2 0
54 1684.5 0
55 1414.7 0
56 1674.5 0
57 1598.7 0
58 1739.1 0
59 1674.6 0
60 1671.8 0
61 1802.0 0
62 1526.8 0
63 1580.9 0
64 1634.8 0
65 1610.3 0
66 1712.0 0
67 1678.8 0
68 1708.1 0
69 1680.6 0
70 2056.0 1
71 1624.0 1
72 2021.4 1
73 1861.1 1
74 1750.8 1
75 1767.5 1
76 1710.3 1
77 2151.5 1
78 2047.9 1
79 1915.4 1
80 1984.7 1
81 1896.5 1
82 2170.8 1
83 2139.9 1
84 2330.5 1
85 2121.8 1
86 2226.8 1
87 1857.9 1
88 2155.9 1
89 2341.7 1
90 2290.2 1
91 2006.5 1
92 2111.9 1
93 1731.3 1
94 1762.2 1
95 1863.2 1
96 1943.5 1
97 1975.2 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D
1589.9 403.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-369.443 -105.751 -1.151 94.649 360.449
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1589.85 18.89 84.16 <2e-16 ***
D 403.59 35.16 11.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 156.9 on 95 degrees of freedom
Multiple R-squared: 0.5811, Adjusted R-squared: 0.5766
F-statistic: 131.8 on 1 and 95 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.654780171 0.69043966 0.3452198
[2,] 0.845806918 0.30838616 0.1541931
[3,] 0.756181904 0.48763619 0.2438181
[4,] 0.681068701 0.63786260 0.3189313
[5,] 0.629655145 0.74068971 0.3703449
[6,] 0.595724996 0.80855001 0.4042750
[7,] 0.502932415 0.99413517 0.4970676
[8,] 0.473012969 0.94602594 0.5269870
[9,] 0.407677817 0.81535563 0.5923222
[10,] 0.343373386 0.68674677 0.6566266
[11,] 0.344786455 0.68957291 0.6552135
[12,] 0.631516219 0.73696756 0.3684838
[13,] 0.556992467 0.88601507 0.4430075
[14,] 0.484710623 0.96942125 0.5152894
[15,] 0.685563167 0.62887367 0.3144368
[16,] 0.659599520 0.68080096 0.3404005
[17,] 0.616286932 0.76742614 0.3837131
[18,] 0.579231459 0.84153708 0.4207685
[19,] 0.514558757 0.97088249 0.4854412
[20,] 0.462892026 0.92578405 0.5371080
[21,] 0.402020125 0.80404025 0.5979799
[22,] 0.341386791 0.68277358 0.6586132
[23,] 0.326353089 0.65270618 0.6736469
[24,] 0.271216083 0.54243217 0.7287839
[25,] 0.224978451 0.44995690 0.7750215
[26,] 0.181124050 0.36224810 0.8188760
[27,] 0.174753693 0.34950739 0.8252463
[28,] 0.137148716 0.27429743 0.8628513
[29,] 0.162809076 0.32561815 0.8371909
[30,] 0.151150527 0.30230105 0.8488495
[31,] 0.126783984 0.25356797 0.8732160
[32,] 0.112162997 0.22432599 0.8878370
[33,] 0.095814017 0.19162803 0.9041860
[34,] 0.081554024 0.16310805 0.9184460
[35,] 0.159293150 0.31858630 0.8407069
[36,] 0.193721293 0.38744259 0.8062787
[37,] 0.164744798 0.32948960 0.8352552
[38,] 0.154925454 0.30985091 0.8450745
[39,] 0.197285804 0.39457161 0.8027142
[40,] 0.223765619 0.44753124 0.7762344
[41,] 0.294551607 0.58910321 0.7054484
[42,] 0.421661800 0.84332360 0.5783382
[43,] 0.377918663 0.75583733 0.6220813
[44,] 0.330518046 0.66103609 0.6694820
[45,] 0.301495423 0.60299085 0.6985046
[46,] 0.266105107 0.53221021 0.7338949
[47,] 0.254626423 0.50925285 0.7453736
[48,] 0.237787768 0.47557554 0.7622122
[49,] 0.198431904 0.39686381 0.8015681
[50,] 0.170115995 0.34023199 0.8298840
[51,] 0.196459926 0.39291985 0.8035401
[52,] 0.165570297 0.33114059 0.8344297
[53,] 0.134493933 0.26898787 0.8655061
[54,] 0.123052616 0.24610523 0.8769474
[55,] 0.099884315 0.19976863 0.9001157
[56,] 0.079602262 0.15920452 0.9203977
[57,] 0.087930203 0.17586041 0.9120698
[58,] 0.073583442 0.14716688 0.9264166
[59,] 0.057047947 0.11409589 0.9429521
[60,] 0.042632101 0.08526420 0.9573679
[61,] 0.031830285 0.06366057 0.9681697
[62,] 0.024708384 0.04941677 0.9752916
[63,] 0.017970908 0.03594182 0.9820291
[64,] 0.013362864 0.02672573 0.9866371
[65,] 0.009342481 0.01868496 0.9906575
[66,] 0.006253316 0.01250663 0.9937467
[67,] 0.024949874 0.04989975 0.9750501
[68,] 0.019095151 0.03819030 0.9809048
[69,] 0.015479670 0.03095934 0.9845203
[70,] 0.021449723 0.04289945 0.9785503
[71,] 0.028400225 0.05680045 0.9715998
[72,] 0.059864523 0.11972905 0.9401355
[73,] 0.066961226 0.13392245 0.9330388
[74,] 0.052019434 0.10403887 0.9479806
[75,] 0.042204584 0.08440917 0.9577954
[76,] 0.030312644 0.06062529 0.9696874
[77,] 0.026319856 0.05263971 0.9736801
[78,] 0.025214324 0.05042865 0.9747857
[79,] 0.020294566 0.04058913 0.9797054
[80,] 0.054404405 0.10880881 0.9455956
[81,] 0.040995138 0.08199028 0.9590049
[82,] 0.051591728 0.10318346 0.9484083
[83,] 0.043543966 0.08708793 0.9564560
[84,] 0.035988119 0.07197624 0.9640119
[85,] 0.147562262 0.29512452 0.8524377
[86,] 0.482387390 0.96477478 0.5176126
[87,] 0.391819835 0.78363967 0.6081802
[88,] 0.592681663 0.81463667 0.4073183
> postscript(file="/var/www/html/freestat/rcomp/tmp/1uzwe1227192425.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/freestat/rcomp/tmp/2id5z1227192425.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/freestat/rcomp/tmp/3rdio1227192425.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/freestat/rcomp/tmp/403q71227192425.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/freestat/rcomp/tmp/57uso1227192425.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 = 97
Frequency = 1
1 2 3 4 5 6
3.149275 -111.950725 143.849275 -20.150725 253.849275 360.449275
7 8 9 10 11 12
67.649275 182.249275 -21.550725 219.949275 56.849275 218.649275
13 14 15 16 17 18
174.049275 35.649275 -51.050725 -247.450725 55.249275 30.049275
19 20 21 22 23 24
-251.750725 -84.350725 -60.750725 -77.950725 66.849275 104.549275
25 26 27 28 29 30
72.449275 -1.150725 -106.550725 -4.250725 69.049275 -5.450725
31 32 33 34 35 36
-119.250725 28.849275 -182.250725 -115.950725 -74.550725 -104.450725
37 38 39 40 41 42
-93.750725 -96.350725 -291.450725 -214.550725 -81.950725 -134.550725
43 44 45 46 47 48
-226.550725 -197.050725 -241.050725 290.449275 79.349275 -46.250725
49 50 51 52 53 54
111.349275 -73.350725 -123.050725 -105.750725 -12.650725 94.649275
55 56 57 58 59 60
-175.150725 84.649275 8.849275 149.249275 84.749275 81.949275
61 62 63 64 65 66
212.149275 -63.050725 -8.950725 44.949275 20.449275 122.149275
67 68 69 70 71 72
88.949275 118.249275 90.749275 62.557143 -369.442857 27.957143
73 74 75 76 77 78
-132.342857 -242.642857 -225.942857 -283.142857 158.057143 54.457143
79 80 81 82 83 84
-78.042857 -8.742857 -96.942857 177.357143 146.457143 337.057143
85 86 87 88 89 90
128.357143 233.357143 -135.542857 162.457143 348.257143 296.757143
91 92 93 94 95 96
13.057143 118.457143 -262.142857 -231.242857 -130.242857 -49.942857
97
-18.242857
> postscript(file="/var/www/html/freestat/rcomp/tmp/6t77i1227192425.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 = 97
Frequency = 1
lag(myerror, k = 1) myerror
0 3.149275 NA
1 -111.950725 3.149275
2 143.849275 -111.950725
3 -20.150725 143.849275
4 253.849275 -20.150725
5 360.449275 253.849275
6 67.649275 360.449275
7 182.249275 67.649275
8 -21.550725 182.249275
9 219.949275 -21.550725
10 56.849275 219.949275
11 218.649275 56.849275
12 174.049275 218.649275
13 35.649275 174.049275
14 -51.050725 35.649275
15 -247.450725 -51.050725
16 55.249275 -247.450725
17 30.049275 55.249275
18 -251.750725 30.049275
19 -84.350725 -251.750725
20 -60.750725 -84.350725
21 -77.950725 -60.750725
22 66.849275 -77.950725
23 104.549275 66.849275
24 72.449275 104.549275
25 -1.150725 72.449275
26 -106.550725 -1.150725
27 -4.250725 -106.550725
28 69.049275 -4.250725
29 -5.450725 69.049275
30 -119.250725 -5.450725
31 28.849275 -119.250725
32 -182.250725 28.849275
33 -115.950725 -182.250725
34 -74.550725 -115.950725
35 -104.450725 -74.550725
36 -93.750725 -104.450725
37 -96.350725 -93.750725
38 -291.450725 -96.350725
39 -214.550725 -291.450725
40 -81.950725 -214.550725
41 -134.550725 -81.950725
42 -226.550725 -134.550725
43 -197.050725 -226.550725
44 -241.050725 -197.050725
45 290.449275 -241.050725
46 79.349275 290.449275
47 -46.250725 79.349275
48 111.349275 -46.250725
49 -73.350725 111.349275
50 -123.050725 -73.350725
51 -105.750725 -123.050725
52 -12.650725 -105.750725
53 94.649275 -12.650725
54 -175.150725 94.649275
55 84.649275 -175.150725
56 8.849275 84.649275
57 149.249275 8.849275
58 84.749275 149.249275
59 81.949275 84.749275
60 212.149275 81.949275
61 -63.050725 212.149275
62 -8.950725 -63.050725
63 44.949275 -8.950725
64 20.449275 44.949275
65 122.149275 20.449275
66 88.949275 122.149275
67 118.249275 88.949275
68 90.749275 118.249275
69 62.557143 90.749275
70 -369.442857 62.557143
71 27.957143 -369.442857
72 -132.342857 27.957143
73 -242.642857 -132.342857
74 -225.942857 -242.642857
75 -283.142857 -225.942857
76 158.057143 -283.142857
77 54.457143 158.057143
78 -78.042857 54.457143
79 -8.742857 -78.042857
80 -96.942857 -8.742857
81 177.357143 -96.942857
82 146.457143 177.357143
83 337.057143 146.457143
84 128.357143 337.057143
85 233.357143 128.357143
86 -135.542857 233.357143
87 162.457143 -135.542857
88 348.257143 162.457143
89 296.757143 348.257143
90 13.057143 296.757143
91 118.457143 13.057143
92 -262.142857 118.457143
93 -231.242857 -262.142857
94 -130.242857 -231.242857
95 -49.942857 -130.242857
96 -18.242857 -49.942857
97 NA -18.242857
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -111.950725 3.149275
[2,] 143.849275 -111.950725
[3,] -20.150725 143.849275
[4,] 253.849275 -20.150725
[5,] 360.449275 253.849275
[6,] 67.649275 360.449275
[7,] 182.249275 67.649275
[8,] -21.550725 182.249275
[9,] 219.949275 -21.550725
[10,] 56.849275 219.949275
[11,] 218.649275 56.849275
[12,] 174.049275 218.649275
[13,] 35.649275 174.049275
[14,] -51.050725 35.649275
[15,] -247.450725 -51.050725
[16,] 55.249275 -247.450725
[17,] 30.049275 55.249275
[18,] -251.750725 30.049275
[19,] -84.350725 -251.750725
[20,] -60.750725 -84.350725
[21,] -77.950725 -60.750725
[22,] 66.849275 -77.950725
[23,] 104.549275 66.849275
[24,] 72.449275 104.549275
[25,] -1.150725 72.449275
[26,] -106.550725 -1.150725
[27,] -4.250725 -106.550725
[28,] 69.049275 -4.250725
[29,] -5.450725 69.049275
[30,] -119.250725 -5.450725
[31,] 28.849275 -119.250725
[32,] -182.250725 28.849275
[33,] -115.950725 -182.250725
[34,] -74.550725 -115.950725
[35,] -104.450725 -74.550725
[36,] -93.750725 -104.450725
[37,] -96.350725 -93.750725
[38,] -291.450725 -96.350725
[39,] -214.550725 -291.450725
[40,] -81.950725 -214.550725
[41,] -134.550725 -81.950725
[42,] -226.550725 -134.550725
[43,] -197.050725 -226.550725
[44,] -241.050725 -197.050725
[45,] 290.449275 -241.050725
[46,] 79.349275 290.449275
[47,] -46.250725 79.349275
[48,] 111.349275 -46.250725
[49,] -73.350725 111.349275
[50,] -123.050725 -73.350725
[51,] -105.750725 -123.050725
[52,] -12.650725 -105.750725
[53,] 94.649275 -12.650725
[54,] -175.150725 94.649275
[55,] 84.649275 -175.150725
[56,] 8.849275 84.649275
[57,] 149.249275 8.849275
[58,] 84.749275 149.249275
[59,] 81.949275 84.749275
[60,] 212.149275 81.949275
[61,] -63.050725 212.149275
[62,] -8.950725 -63.050725
[63,] 44.949275 -8.950725
[64,] 20.449275 44.949275
[65,] 122.149275 20.449275
[66,] 88.949275 122.149275
[67,] 118.249275 88.949275
[68,] 90.749275 118.249275
[69,] 62.557143 90.749275
[70,] -369.442857 62.557143
[71,] 27.957143 -369.442857
[72,] -132.342857 27.957143
[73,] -242.642857 -132.342857
[74,] -225.942857 -242.642857
[75,] -283.142857 -225.942857
[76,] 158.057143 -283.142857
[77,] 54.457143 158.057143
[78,] -78.042857 54.457143
[79,] -8.742857 -78.042857
[80,] -96.942857 -8.742857
[81,] 177.357143 -96.942857
[82,] 146.457143 177.357143
[83,] 337.057143 146.457143
[84,] 128.357143 337.057143
[85,] 233.357143 128.357143
[86,] -135.542857 233.357143
[87,] 162.457143 -135.542857
[88,] 348.257143 162.457143
[89,] 296.757143 348.257143
[90,] 13.057143 296.757143
[91,] 118.457143 13.057143
[92,] -262.142857 118.457143
[93,] -231.242857 -262.142857
[94,] -130.242857 -231.242857
[95,] -49.942857 -130.242857
[96,] -18.242857 -49.942857
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -111.950725 3.149275
2 143.849275 -111.950725
3 -20.150725 143.849275
4 253.849275 -20.150725
5 360.449275 253.849275
6 67.649275 360.449275
7 182.249275 67.649275
8 -21.550725 182.249275
9 219.949275 -21.550725
10 56.849275 219.949275
11 218.649275 56.849275
12 174.049275 218.649275
13 35.649275 174.049275
14 -51.050725 35.649275
15 -247.450725 -51.050725
16 55.249275 -247.450725
17 30.049275 55.249275
18 -251.750725 30.049275
19 -84.350725 -251.750725
20 -60.750725 -84.350725
21 -77.950725 -60.750725
22 66.849275 -77.950725
23 104.549275 66.849275
24 72.449275 104.549275
25 -1.150725 72.449275
26 -106.550725 -1.150725
27 -4.250725 -106.550725
28 69.049275 -4.250725
29 -5.450725 69.049275
30 -119.250725 -5.450725
31 28.849275 -119.250725
32 -182.250725 28.849275
33 -115.950725 -182.250725
34 -74.550725 -115.950725
35 -104.450725 -74.550725
36 -93.750725 -104.450725
37 -96.350725 -93.750725
38 -291.450725 -96.350725
39 -214.550725 -291.450725
40 -81.950725 -214.550725
41 -134.550725 -81.950725
42 -226.550725 -134.550725
43 -197.050725 -226.550725
44 -241.050725 -197.050725
45 290.449275 -241.050725
46 79.349275 290.449275
47 -46.250725 79.349275
48 111.349275 -46.250725
49 -73.350725 111.349275
50 -123.050725 -73.350725
51 -105.750725 -123.050725
52 -12.650725 -105.750725
53 94.649275 -12.650725
54 -175.150725 94.649275
55 84.649275 -175.150725
56 8.849275 84.649275
57 149.249275 8.849275
58 84.749275 149.249275
59 81.949275 84.749275
60 212.149275 81.949275
61 -63.050725 212.149275
62 -8.950725 -63.050725
63 44.949275 -8.950725
64 20.449275 44.949275
65 122.149275 20.449275
66 88.949275 122.149275
67 118.249275 88.949275
68 90.749275 118.249275
69 62.557143 90.749275
70 -369.442857 62.557143
71 27.957143 -369.442857
72 -132.342857 27.957143
73 -242.642857 -132.342857
74 -225.942857 -242.642857
75 -283.142857 -225.942857
76 158.057143 -283.142857
77 54.457143 158.057143
78 -78.042857 54.457143
79 -8.742857 -78.042857
80 -96.942857 -8.742857
81 177.357143 -96.942857
82 146.457143 177.357143
83 337.057143 146.457143
84 128.357143 337.057143
85 233.357143 128.357143
86 -135.542857 233.357143
87 162.457143 -135.542857
88 348.257143 162.457143
89 296.757143 348.257143
90 13.057143 296.757143
91 118.457143 13.057143
92 -262.142857 118.457143
93 -231.242857 -262.142857
94 -130.242857 -231.242857
95 -49.942857 -130.242857
96 -18.242857 -49.942857
> 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/freestat/rcomp/tmp/7saz01227192425.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/freestat/rcomp/tmp/8ue8w1227192425.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/freestat/rcomp/tmp/9bdgv1227192425.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/freestat/rcomp/tmp/10mf561227192425.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11er2b1227192425.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/freestat/rcomp/tmp/12ysh01227192425.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/freestat/rcomp/tmp/13daqm1227192425.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/freestat/rcomp/tmp/14ei5i1227192425.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/freestat/rcomp/tmp/15a4ql1227192425.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/freestat/rcomp/tmp/16ox0v1227192425.tab")
+ }
>
> system("convert tmp/1uzwe1227192425.ps tmp/1uzwe1227192425.png")
> system("convert tmp/2id5z1227192425.ps tmp/2id5z1227192425.png")
> system("convert tmp/3rdio1227192425.ps tmp/3rdio1227192425.png")
> system("convert tmp/403q71227192425.ps tmp/403q71227192425.png")
> system("convert tmp/57uso1227192425.ps tmp/57uso1227192425.png")
> system("convert tmp/6t77i1227192425.ps tmp/6t77i1227192425.png")
> system("convert tmp/7saz01227192425.ps tmp/7saz01227192425.png")
> system("convert tmp/8ue8w1227192425.ps tmp/8ue8w1227192425.png")
> system("convert tmp/9bdgv1227192425.ps tmp/9bdgv1227192425.png")
> system("convert tmp/10mf561227192425.ps tmp/10mf561227192425.png")
>
>
> proc.time()
user system elapsed
4.506 2.698 6.365