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(101.02,0,100.67,0,100.47,0,100.38,0,100.33,0,100.34,0,100.37,0,100.39,0,100.21,0,100.21,0,100.22,0,100.28,0,100.25,0,100.25,0,100.21,0,100.16,0,100.18,0,100.1,1,99.96,1,99.88,1,99.88,1,99.86,1,99.84,1,99.8,1,99.82,1,99.81,1,99.92,1,100.03,1,99.99,1,100.02,1,100.01,1,100.13,1,100.33,1,100.13,1,99.96,1,100.05,1,99.83,1,99.8,1,100.01,1,100.1,1,100.13,1,100.16,1,100.41,1,101.34,1,101.65,1,101.85,1,102.07,1,102.12,1,102.14,1,102.21,1,102.28,1,102.19,1,102.33,1,102.54,1,102.44,1,102.78,1,102.9,1,103.08,1,102.77,1,102.65,1,102.71,1,103.29,1,102.86,1,103.45,1,103.72,1,103.65,1,103.83,1,104.45,1,105.14,1,105.07,1,105.31,1,105.19,1,105.3,1,105.02,1,105.17,1,105.28,1,105.45,1,105.38,1,105.8,1,105.96,1,105.08,1,105.11,1,105.61,1,105.5,1),dim=c(2,84),dimnames=list(c('Suiker','Dummy'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Suiker','Dummy'),1:84))
> 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
Suiker Dummy
1 101.02 0
2 100.67 0
3 100.47 0
4 100.38 0
5 100.33 0
6 100.34 0
7 100.37 0
8 100.39 0
9 100.21 0
10 100.21 0
11 100.22 0
12 100.28 0
13 100.25 0
14 100.25 0
15 100.21 0
16 100.16 0
17 100.18 0
18 100.10 1
19 99.96 1
20 99.88 1
21 99.88 1
22 99.86 1
23 99.84 1
24 99.80 1
25 99.82 1
26 99.81 1
27 99.92 1
28 100.03 1
29 99.99 1
30 100.02 1
31 100.01 1
32 100.13 1
33 100.33 1
34 100.13 1
35 99.96 1
36 100.05 1
37 99.83 1
38 99.80 1
39 100.01 1
40 100.10 1
41 100.13 1
42 100.16 1
43 100.41 1
44 101.34 1
45 101.65 1
46 101.85 1
47 102.07 1
48 102.12 1
49 102.14 1
50 102.21 1
51 102.28 1
52 102.19 1
53 102.33 1
54 102.54 1
55 102.44 1
56 102.78 1
57 102.90 1
58 103.08 1
59 102.77 1
60 102.65 1
61 102.71 1
62 103.29 1
63 102.86 1
64 103.45 1
65 103.72 1
66 103.65 1
67 103.83 1
68 104.45 1
69 105.14 1
70 105.07 1
71 105.31 1
72 105.19 1
73 105.30 1
74 105.02 1
75 105.17 1
76 105.28 1
77 105.45 1
78 105.38 1
79 105.80 1
80 105.96 1
81 105.08 1
82 105.11 1
83 105.61 1
84 105.50 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy
100.349 1.929
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.47806 -2.14806 -0.06874 1.05194 3.68194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.3494 0.4674 214.713 < 2e-16 ***
Dummy 1.9286 0.5233 3.685 0.000408 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.927 on 82 degrees of freedom
Multiple R-squared: 0.1421, Adjusted R-squared: 0.1316
F-statistic: 13.58 on 1 and 82 DF, p-value: 0.0004085
> 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,] 5.940067e-03 1.188013e-02 9.940599e-01
[2,] 1.017828e-03 2.035655e-03 9.989822e-01
[3,] 1.501725e-04 3.003449e-04 9.998498e-01
[4,] 1.976261e-05 3.952521e-05 9.999802e-01
[5,] 3.965222e-06 7.930444e-06 9.999960e-01
[6,] 7.160360e-07 1.432072e-06 9.999993e-01
[7,] 1.154890e-07 2.309780e-07 9.999999e-01
[8,] 1.497264e-08 2.994528e-08 1.000000e+00
[9,] 1.978207e-09 3.956414e-09 1.000000e+00
[10,] 2.486807e-10 4.973614e-10 1.000000e+00
[11,] 3.342725e-11 6.685450e-11 1.000000e+00
[12,] 5.104985e-12 1.020997e-11 1.000000e+00
[13,] 6.860170e-13 1.372034e-12 1.000000e+00
[14,] 7.504573e-14 1.500915e-13 1.000000e+00
[15,] 9.211061e-15 1.842212e-14 1.000000e+00
[16,] 1.220882e-15 2.441764e-15 1.000000e+00
[17,] 1.497361e-16 2.994723e-16 1.000000e+00
[18,] 1.856098e-17 3.712196e-17 1.000000e+00
[19,] 2.359430e-18 4.718860e-18 1.000000e+00
[20,] 3.294430e-19 6.588859e-19 1.000000e+00
[21,] 4.317311e-20 8.634622e-20 1.000000e+00
[22,] 5.851566e-21 1.170313e-20 1.000000e+00
[23,] 7.399025e-22 1.479805e-21 1.000000e+00
[24,] 1.228292e-22 2.456584e-22 1.000000e+00
[25,] 1.815038e-23 3.630077e-23 1.000000e+00
[26,] 2.990206e-24 5.980411e-24 1.000000e+00
[27,] 4.951187e-25 9.902374e-25 1.000000e+00
[28,] 1.427352e-25 2.854704e-25 1.000000e+00
[29,] 2.073253e-25 4.146507e-25 1.000000e+00
[30,] 5.647495e-26 1.129499e-25 1.000000e+00
[31,] 1.255368e-26 2.510736e-26 1.000000e+00
[32,] 3.270422e-27 6.540845e-27 1.000000e+00
[33,] 1.447999e-27 2.895997e-27 1.000000e+00
[34,] 9.193389e-28 1.838678e-27 1.000000e+00
[35,] 4.275443e-28 8.550887e-28 1.000000e+00
[36,] 3.120464e-28 6.240928e-28 1.000000e+00
[37,] 3.475892e-28 6.951785e-28 1.000000e+00
[38,] 6.517463e-28 1.303493e-27 1.000000e+00
[39,] 1.264666e-26 2.529332e-26 1.000000e+00
[40,] 8.215665e-20 1.643133e-19 1.000000e+00
[41,] 4.554788e-15 9.109575e-15 1.000000e+00
[42,] 8.489045e-12 1.697809e-11 1.000000e+00
[43,] 2.587276e-09 5.174551e-09 1.000000e+00
[44,] 1.355701e-07 2.711402e-07 9.999999e-01
[45,] 2.491001e-06 4.982002e-06 9.999975e-01
[46,] 2.585099e-05 5.170198e-05 9.999741e-01
[47,] 1.770746e-04 3.541493e-04 9.998229e-01
[48,] 8.253400e-04 1.650680e-03 9.991747e-01
[49,] 3.273636e-03 6.547272e-03 9.967264e-01
[50,] 1.095743e-02 2.191485e-02 9.890426e-01
[51,] 3.009934e-02 6.019868e-02 9.699007e-01
[52,] 7.019848e-02 1.403970e-01 9.298015e-01
[53,] 1.372545e-01 2.745089e-01 8.627455e-01
[54,] 2.295158e-01 4.590316e-01 7.704842e-01
[55,] 3.612592e-01 7.225183e-01 6.387408e-01
[56,] 5.467285e-01 9.065431e-01 4.532715e-01
[57,] 7.490683e-01 5.018634e-01 2.509317e-01
[58,] 8.607271e-01 2.785457e-01 1.392729e-01
[59,] 9.662138e-01 6.757238e-02 3.378619e-02
[60,] 9.912818e-01 1.743632e-02 8.718161e-03
[61,] 9.979316e-01 4.136709e-03 2.068354e-03
[62,] 9.998480e-01 3.039282e-04 1.519641e-04
[63,] 9.999980e-01 3.967361e-06 1.983680e-06
[64,] 9.999999e-01 2.649891e-07 1.324945e-07
[65,] 9.999998e-01 4.492206e-07 2.246103e-07
[66,] 9.999996e-01 7.498937e-07 3.749468e-07
[67,] 9.999989e-01 2.168568e-06 1.084284e-06
[68,] 9.999969e-01 6.132768e-06 3.066384e-06
[69,] 9.999893e-01 2.144978e-05 1.072489e-05
[70,] 9.999788e-01 4.238136e-05 2.119068e-05
[71,] 9.999350e-01 1.300205e-04 6.501026e-05
[72,] 9.997456e-01 5.087861e-04 2.543931e-04
[73,] 9.988263e-01 2.347498e-03 1.173749e-03
[74,] 9.948628e-01 1.027433e-02 5.137167e-03
[75,] 9.826201e-01 3.475979e-02 1.737990e-02
> postscript(file="/var/www/html/rcomp/tmp/1a8tk1229364603.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/2qa1g1229364603.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/3pxt71229364603.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/498gd1229364603.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/5kvvu1229364603.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 = 84
Frequency = 1
1 2 3 4 5 6
0.670588235 0.320588235 0.120588235 0.030588235 -0.019411765 -0.009411765
7 8 9 10 11 12
0.020588235 0.040588235 -0.139411765 -0.139411765 -0.129411765 -0.069411765
13 14 15 16 17 18
-0.099411765 -0.099411765 -0.139411765 -0.189411765 -0.169411765 -2.178059701
19 20 21 22 23 24
-2.318059701 -2.398059701 -2.398059701 -2.418059701 -2.438059701 -2.478059701
25 26 27 28 29 30
-2.458059701 -2.468059701 -2.358059701 -2.248059701 -2.288059701 -2.258059701
31 32 33 34 35 36
-2.268059701 -2.148059701 -1.948059701 -2.148059701 -2.318059701 -2.228059701
37 38 39 40 41 42
-2.448059701 -2.478059701 -2.268059701 -2.178059701 -2.148059701 -2.118059701
43 44 45 46 47 48
-1.868059701 -0.938059701 -0.628059701 -0.428059701 -0.208059701 -0.158059701
49 50 51 52 53 54
-0.138059701 -0.068059701 0.001940299 -0.088059701 0.051940299 0.261940299
55 56 57 58 59 60
0.161940299 0.501940299 0.621940299 0.801940299 0.491940299 0.371940299
61 62 63 64 65 66
0.431940299 1.011940299 0.581940299 1.171940299 1.441940299 1.371940299
67 68 69 70 71 72
1.551940299 2.171940299 2.861940299 2.791940299 3.031940299 2.911940299
73 74 75 76 77 78
3.021940299 2.741940299 2.891940299 3.001940299 3.171940299 3.101940299
79 80 81 82 83 84
3.521940299 3.681940299 2.801940299 2.831940299 3.331940299 3.221940299
> postscript(file="/var/www/html/rcomp/tmp/6qd041229364603.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 0.670588235 NA
1 0.320588235 0.670588235
2 0.120588235 0.320588235
3 0.030588235 0.120588235
4 -0.019411765 0.030588235
5 -0.009411765 -0.019411765
6 0.020588235 -0.009411765
7 0.040588235 0.020588235
8 -0.139411765 0.040588235
9 -0.139411765 -0.139411765
10 -0.129411765 -0.139411765
11 -0.069411765 -0.129411765
12 -0.099411765 -0.069411765
13 -0.099411765 -0.099411765
14 -0.139411765 -0.099411765
15 -0.189411765 -0.139411765
16 -0.169411765 -0.189411765
17 -2.178059701 -0.169411765
18 -2.318059701 -2.178059701
19 -2.398059701 -2.318059701
20 -2.398059701 -2.398059701
21 -2.418059701 -2.398059701
22 -2.438059701 -2.418059701
23 -2.478059701 -2.438059701
24 -2.458059701 -2.478059701
25 -2.468059701 -2.458059701
26 -2.358059701 -2.468059701
27 -2.248059701 -2.358059701
28 -2.288059701 -2.248059701
29 -2.258059701 -2.288059701
30 -2.268059701 -2.258059701
31 -2.148059701 -2.268059701
32 -1.948059701 -2.148059701
33 -2.148059701 -1.948059701
34 -2.318059701 -2.148059701
35 -2.228059701 -2.318059701
36 -2.448059701 -2.228059701
37 -2.478059701 -2.448059701
38 -2.268059701 -2.478059701
39 -2.178059701 -2.268059701
40 -2.148059701 -2.178059701
41 -2.118059701 -2.148059701
42 -1.868059701 -2.118059701
43 -0.938059701 -1.868059701
44 -0.628059701 -0.938059701
45 -0.428059701 -0.628059701
46 -0.208059701 -0.428059701
47 -0.158059701 -0.208059701
48 -0.138059701 -0.158059701
49 -0.068059701 -0.138059701
50 0.001940299 -0.068059701
51 -0.088059701 0.001940299
52 0.051940299 -0.088059701
53 0.261940299 0.051940299
54 0.161940299 0.261940299
55 0.501940299 0.161940299
56 0.621940299 0.501940299
57 0.801940299 0.621940299
58 0.491940299 0.801940299
59 0.371940299 0.491940299
60 0.431940299 0.371940299
61 1.011940299 0.431940299
62 0.581940299 1.011940299
63 1.171940299 0.581940299
64 1.441940299 1.171940299
65 1.371940299 1.441940299
66 1.551940299 1.371940299
67 2.171940299 1.551940299
68 2.861940299 2.171940299
69 2.791940299 2.861940299
70 3.031940299 2.791940299
71 2.911940299 3.031940299
72 3.021940299 2.911940299
73 2.741940299 3.021940299
74 2.891940299 2.741940299
75 3.001940299 2.891940299
76 3.171940299 3.001940299
77 3.101940299 3.171940299
78 3.521940299 3.101940299
79 3.681940299 3.521940299
80 2.801940299 3.681940299
81 2.831940299 2.801940299
82 3.331940299 2.831940299
83 3.221940299 3.331940299
84 NA 3.221940299
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.320588235 0.670588235
[2,] 0.120588235 0.320588235
[3,] 0.030588235 0.120588235
[4,] -0.019411765 0.030588235
[5,] -0.009411765 -0.019411765
[6,] 0.020588235 -0.009411765
[7,] 0.040588235 0.020588235
[8,] -0.139411765 0.040588235
[9,] -0.139411765 -0.139411765
[10,] -0.129411765 -0.139411765
[11,] -0.069411765 -0.129411765
[12,] -0.099411765 -0.069411765
[13,] -0.099411765 -0.099411765
[14,] -0.139411765 -0.099411765
[15,] -0.189411765 -0.139411765
[16,] -0.169411765 -0.189411765
[17,] -2.178059701 -0.169411765
[18,] -2.318059701 -2.178059701
[19,] -2.398059701 -2.318059701
[20,] -2.398059701 -2.398059701
[21,] -2.418059701 -2.398059701
[22,] -2.438059701 -2.418059701
[23,] -2.478059701 -2.438059701
[24,] -2.458059701 -2.478059701
[25,] -2.468059701 -2.458059701
[26,] -2.358059701 -2.468059701
[27,] -2.248059701 -2.358059701
[28,] -2.288059701 -2.248059701
[29,] -2.258059701 -2.288059701
[30,] -2.268059701 -2.258059701
[31,] -2.148059701 -2.268059701
[32,] -1.948059701 -2.148059701
[33,] -2.148059701 -1.948059701
[34,] -2.318059701 -2.148059701
[35,] -2.228059701 -2.318059701
[36,] -2.448059701 -2.228059701
[37,] -2.478059701 -2.448059701
[38,] -2.268059701 -2.478059701
[39,] -2.178059701 -2.268059701
[40,] -2.148059701 -2.178059701
[41,] -2.118059701 -2.148059701
[42,] -1.868059701 -2.118059701
[43,] -0.938059701 -1.868059701
[44,] -0.628059701 -0.938059701
[45,] -0.428059701 -0.628059701
[46,] -0.208059701 -0.428059701
[47,] -0.158059701 -0.208059701
[48,] -0.138059701 -0.158059701
[49,] -0.068059701 -0.138059701
[50,] 0.001940299 -0.068059701
[51,] -0.088059701 0.001940299
[52,] 0.051940299 -0.088059701
[53,] 0.261940299 0.051940299
[54,] 0.161940299 0.261940299
[55,] 0.501940299 0.161940299
[56,] 0.621940299 0.501940299
[57,] 0.801940299 0.621940299
[58,] 0.491940299 0.801940299
[59,] 0.371940299 0.491940299
[60,] 0.431940299 0.371940299
[61,] 1.011940299 0.431940299
[62,] 0.581940299 1.011940299
[63,] 1.171940299 0.581940299
[64,] 1.441940299 1.171940299
[65,] 1.371940299 1.441940299
[66,] 1.551940299 1.371940299
[67,] 2.171940299 1.551940299
[68,] 2.861940299 2.171940299
[69,] 2.791940299 2.861940299
[70,] 3.031940299 2.791940299
[71,] 2.911940299 3.031940299
[72,] 3.021940299 2.911940299
[73,] 2.741940299 3.021940299
[74,] 2.891940299 2.741940299
[75,] 3.001940299 2.891940299
[76,] 3.171940299 3.001940299
[77,] 3.101940299 3.171940299
[78,] 3.521940299 3.101940299
[79,] 3.681940299 3.521940299
[80,] 2.801940299 3.681940299
[81,] 2.831940299 2.801940299
[82,] 3.331940299 2.831940299
[83,] 3.221940299 3.331940299
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.320588235 0.670588235
2 0.120588235 0.320588235
3 0.030588235 0.120588235
4 -0.019411765 0.030588235
5 -0.009411765 -0.019411765
6 0.020588235 -0.009411765
7 0.040588235 0.020588235
8 -0.139411765 0.040588235
9 -0.139411765 -0.139411765
10 -0.129411765 -0.139411765
11 -0.069411765 -0.129411765
12 -0.099411765 -0.069411765
13 -0.099411765 -0.099411765
14 -0.139411765 -0.099411765
15 -0.189411765 -0.139411765
16 -0.169411765 -0.189411765
17 -2.178059701 -0.169411765
18 -2.318059701 -2.178059701
19 -2.398059701 -2.318059701
20 -2.398059701 -2.398059701
21 -2.418059701 -2.398059701
22 -2.438059701 -2.418059701
23 -2.478059701 -2.438059701
24 -2.458059701 -2.478059701
25 -2.468059701 -2.458059701
26 -2.358059701 -2.468059701
27 -2.248059701 -2.358059701
28 -2.288059701 -2.248059701
29 -2.258059701 -2.288059701
30 -2.268059701 -2.258059701
31 -2.148059701 -2.268059701
32 -1.948059701 -2.148059701
33 -2.148059701 -1.948059701
34 -2.318059701 -2.148059701
35 -2.228059701 -2.318059701
36 -2.448059701 -2.228059701
37 -2.478059701 -2.448059701
38 -2.268059701 -2.478059701
39 -2.178059701 -2.268059701
40 -2.148059701 -2.178059701
41 -2.118059701 -2.148059701
42 -1.868059701 -2.118059701
43 -0.938059701 -1.868059701
44 -0.628059701 -0.938059701
45 -0.428059701 -0.628059701
46 -0.208059701 -0.428059701
47 -0.158059701 -0.208059701
48 -0.138059701 -0.158059701
49 -0.068059701 -0.138059701
50 0.001940299 -0.068059701
51 -0.088059701 0.001940299
52 0.051940299 -0.088059701
53 0.261940299 0.051940299
54 0.161940299 0.261940299
55 0.501940299 0.161940299
56 0.621940299 0.501940299
57 0.801940299 0.621940299
58 0.491940299 0.801940299
59 0.371940299 0.491940299
60 0.431940299 0.371940299
61 1.011940299 0.431940299
62 0.581940299 1.011940299
63 1.171940299 0.581940299
64 1.441940299 1.171940299
65 1.371940299 1.441940299
66 1.551940299 1.371940299
67 2.171940299 1.551940299
68 2.861940299 2.171940299
69 2.791940299 2.861940299
70 3.031940299 2.791940299
71 2.911940299 3.031940299
72 3.021940299 2.911940299
73 2.741940299 3.021940299
74 2.891940299 2.741940299
75 3.001940299 2.891940299
76 3.171940299 3.001940299
77 3.101940299 3.171940299
78 3.521940299 3.101940299
79 3.681940299 3.521940299
80 2.801940299 3.681940299
81 2.831940299 2.801940299
82 3.331940299 2.831940299
83 3.221940299 3.331940299
> 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/79jov1229364603.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/86ual1229364603.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/9xt0r1229364603.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/10euth1229364603.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/11zhmb1229364603.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/128f171229364603.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/138nxp1229364604.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/14o1lh1229364604.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/155pvq1229364604.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/16khri1229364604.tab")
+ }
>
> system("convert tmp/1a8tk1229364603.ps tmp/1a8tk1229364603.png")
> system("convert tmp/2qa1g1229364603.ps tmp/2qa1g1229364603.png")
> system("convert tmp/3pxt71229364603.ps tmp/3pxt71229364603.png")
> system("convert tmp/498gd1229364603.ps tmp/498gd1229364603.png")
> system("convert tmp/5kvvu1229364603.ps tmp/5kvvu1229364603.png")
> system("convert tmp/6qd041229364603.ps tmp/6qd041229364603.png")
> system("convert tmp/79jov1229364603.ps tmp/79jov1229364603.png")
> system("convert tmp/86ual1229364603.ps tmp/86ual1229364603.png")
> system("convert tmp/9xt0r1229364603.ps tmp/9xt0r1229364603.png")
> system("convert tmp/10euth1229364603.ps tmp/10euth1229364603.png")
>
>
> proc.time()
user system elapsed
2.708 1.570 3.726