R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(1966,1,1966,2,1966,3,1966,4,1966,5,1966,6,1966,7,1966,8,1966,9,1966,10,1966,11,1966,12,1967,1,1967,2,1967,3,1967,4,1967,5,1967,6,1967,7,1967,8,1967,9,1967,10,1967,11,1967,12,1968,1,1968,2,1968,3,1968,4,1968,5,1968,6,1968,7,1968,8,1968,9,1968,10,1968,11,1968,12,1969,1,1969,2,1969,3,1969,4,1969,5,1969,6,1969,7,1969,8,1969,9,1969,10,1969,11,1969,12,1970,1,1970,2,1970,3,1970,4,1970,5,1970,6,1970,7,1970,8,1970,9,1970,10,1970,11,1970,12,1971,1,1971,2,1971,3,1971,4,1971,5,1971,6,1971,7,1971,8,1971,9,1971,10,1971,11,1971,12,1972,1,1972,2,1972,3,1972,4,1972,5,1972,6,1972,7,1972,8,1972,9,1972,10,1972,11,1972,12,1973,1,1973,2,1973,3,1973,4,1973,5,1973,6,1973,7,1973,8,1973,9,1973,10,1973,11,1973,12,1974,1,1974,2,1974,3,1974,4,1974,5,1974,6,1974,7,1974,8,1974,9,1974,10,1974,11,1974,12,1975,1,1975,2,1975,3,1975,4,1975,5,1975,6,1975,7,1975,8,1975,9,1975,10),dim=c(2,118),dimnames=list(c('Jaartal','Maand'),1:118))
> y <- array(NA,dim=c(2,118),dimnames=list(c('Jaartal','Maand'),1:118))
> 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
> 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
Jaartal Maand
1 1966 1
2 1966 2
3 1966 3
4 1966 4
5 1966 5
6 1966 6
7 1966 7
8 1966 8
9 1966 9
10 1966 10
11 1966 11
12 1966 12
13 1967 1
14 1967 2
15 1967 3
16 1967 4
17 1967 5
18 1967 6
19 1967 7
20 1967 8
21 1967 9
22 1967 10
23 1967 11
24 1967 12
25 1968 1
26 1968 2
27 1968 3
28 1968 4
29 1968 5
30 1968 6
31 1968 7
32 1968 8
33 1968 9
34 1968 10
35 1968 11
36 1968 12
37 1969 1
38 1969 2
39 1969 3
40 1969 4
41 1969 5
42 1969 6
43 1969 7
44 1969 8
45 1969 9
46 1969 10
47 1969 11
48 1969 12
49 1970 1
50 1970 2
51 1970 3
52 1970 4
53 1970 5
54 1970 6
55 1970 7
56 1970 8
57 1970 9
58 1970 10
59 1970 11
60 1970 12
61 1971 1
62 1971 2
63 1971 3
64 1971 4
65 1971 5
66 1971 6
67 1971 7
68 1971 8
69 1971 9
70 1971 10
71 1971 11
72 1971 12
73 1972 1
74 1972 2
75 1972 3
76 1972 4
77 1972 5
78 1972 6
79 1972 7
80 1972 8
81 1972 9
82 1972 10
83 1972 11
84 1972 12
85 1973 1
86 1973 2
87 1973 3
88 1973 4
89 1973 5
90 1973 6
91 1973 7
92 1973 8
93 1973 9
94 1973 10
95 1973 11
96 1973 12
97 1974 1
98 1974 2
99 1974 3
100 1974 4
101 1974 5
102 1974 6
103 1974 7
104 1974 8
105 1974 9
106 1974 10
107 1974 11
108 1974 12
109 1975 1
110 1975 2
111 1975 3
112 1975 4
113 1975 5
114 1975 6
115 1975 7
116 1975 8
117 1975 9
118 1975 10
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Maand
1970.63668 -0.03319
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4.6035 -2.4292 -0.2549 2.5210 4.6953
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1970.63668 0.55945 3522.437 <2e-16 ***
Maand -0.03319 0.07696 -0.431 0.667
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.858 on 116 degrees of freedom
Multiple R-squared: 0.001601, Adjusted R-squared: -0.007006
F-statistic: 0.186 on 1 and 116 DF, p-value: 0.6671
> 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,] 1.085450e-41 2.170899e-41 1.000000e+00
[2,] 6.469261e-54 1.293852e-53 1.000000e+00
[3,] 1.780784e-72 3.561567e-72 1.000000e+00
[4,] 1.705437e-78 3.410874e-78 1.000000e+00
[5,] 6.941016e-94 1.388203e-93 1.000000e+00
[6,] 2.095607e-107 4.191214e-107 1.000000e+00
[7,] 1.069740e-125 2.139479e-125 1.000000e+00
[8,] 3.066289e-129 6.132577e-129 1.000000e+00
[9,] 7.531378e-10 1.506276e-09 1.000000e+00
[10,] 2.051767e-09 4.103535e-09 1.000000e+00
[11,] 2.166710e-09 4.333421e-09 1.000000e+00
[12,] 1.795862e-09 3.591725e-09 1.000000e+00
[13,] 1.409261e-09 2.818521e-09 1.000000e+00
[14,] 1.106947e-09 2.213895e-09 1.000000e+00
[15,] 8.703029e-10 1.740606e-09 1.000000e+00
[16,] 6.658928e-10 1.331786e-09 1.000000e+00
[17,] 4.773973e-10 9.547946e-10 1.000000e+00
[18,] 3.106508e-10 6.213015e-10 1.000000e+00
[19,] 1.808870e-10 3.617741e-10 1.000000e+00
[20,] 9.508538e-11 1.901708e-10 1.000000e+00
[21,] 5.914491e-10 1.182898e-09 1.000000e+00
[22,] 1.585929e-09 3.171858e-09 1.000000e+00
[23,] 3.121531e-09 6.243062e-09 1.000000e+00
[24,] 5.459450e-09 1.091890e-08 1.000000e+00
[25,] 9.190180e-09 1.838036e-08 1.000000e+00
[26,] 1.535190e-08 3.070380e-08 1.000000e+00
[27,] 2.553321e-08 5.106642e-08 1.000000e+00
[28,] 4.173745e-08 8.347490e-08 1.000000e+00
[29,] 6.567639e-08 1.313528e-07 9.999999e-01
[30,] 9.751975e-08 1.950395e-07 9.999999e-01
[31,] 1.352215e-07 2.704431e-07 9.999999e-01
[32,] 1.763374e-07 3.526749e-07 9.999998e-01
[33,] 7.355951e-07 1.471190e-06 9.999993e-01
[34,] 2.158836e-06 4.317672e-06 9.999978e-01
[35,] 5.386732e-06 1.077346e-05 9.999946e-01
[36,] 1.240743e-05 2.481485e-05 9.999876e-01
[37,] 2.735363e-05 5.470726e-05 9.999726e-01
[38,] 5.848597e-05 1.169719e-04 9.999415e-01
[39,] 1.213575e-04 2.427150e-04 9.998786e-01
[40,] 2.429405e-04 4.858811e-04 9.997571e-01
[41,] 4.655747e-04 9.311494e-04 9.995344e-01
[42,] 8.498270e-04 1.699654e-03 9.991502e-01
[43,] 1.480823e-03 2.961645e-03 9.985192e-01
[44,] 2.499485e-03 4.998970e-03 9.975005e-01
[45,] 6.139382e-03 1.227876e-02 9.938606e-01
[46,] 1.273957e-02 2.547914e-02 9.872604e-01
[47,] 2.391456e-02 4.782912e-02 9.760854e-01
[48,] 4.180346e-02 8.360693e-02 9.581965e-01
[49,] 6.884647e-02 1.376929e-01 9.311535e-01
[50,] 1.072897e-01 2.145794e-01 8.927103e-01
[51,] 1.585006e-01 3.170012e-01 8.414994e-01
[52,] 2.223636e-01 4.447273e-01 7.776364e-01
[53,] 2.971416e-01 5.942833e-01 7.028584e-01
[54,] 3.800789e-01 7.601578e-01 6.199211e-01
[55,] 4.686894e-01 9.373787e-01 5.313106e-01
[56,] 5.623266e-01 8.753468e-01 4.376734e-01
[57,] 6.706271e-01 6.587459e-01 3.293729e-01
[58,] 7.608259e-01 4.783483e-01 2.391741e-01
[59,] 8.334749e-01 3.330501e-01 1.665251e-01
[60,] 8.890560e-01 2.218879e-01 1.109440e-01
[61,] 9.291272e-01 1.417457e-01 7.087283e-02
[62,] 9.563706e-01 8.725884e-02 4.362942e-02
[63,] 9.739565e-01 5.208699e-02 2.604349e-02
[64,] 9.848464e-01 3.030722e-02 1.515361e-02
[65,] 9.913931e-01 1.721376e-02 8.606879e-03
[66,] 9.952609e-01 9.478194e-03 4.739097e-03
[67,] 9.975294e-01 4.941213e-03 2.470606e-03
[68,] 9.988505e-01 2.299018e-03 1.149509e-03
[69,] 9.994199e-01 1.160269e-03 5.801344e-04
[70,] 9.997143e-01 5.713562e-04 2.856781e-04
[71,] 9.998638e-01 2.724730e-04 1.362365e-04
[72,] 9.999366e-01 1.268146e-04 6.340728e-05
[73,] 9.999707e-01 5.853103e-05 2.926551e-05
[74,] 9.999864e-01 2.726622e-05 1.363311e-05
[75,] 9.999935e-01 1.299391e-05 6.496957e-06
[76,] 9.999968e-01 6.363685e-06 3.181842e-06
[77,] 9.999984e-01 3.173443e-06 1.586722e-06
[78,] 9.999992e-01 1.562612e-06 7.813062e-07
[79,] 9.999996e-01 7.082648e-07 3.541324e-07
[80,] 9.999999e-01 2.519019e-07 1.259509e-07
[81,] 9.999999e-01 1.654059e-07 8.270294e-08
[82,] 9.999999e-01 1.027357e-07 5.136787e-08
[83,] 1.000000e+00 6.071822e-08 3.035911e-08
[84,] 1.000000e+00 3.508328e-08 1.754164e-08
[85,] 1.000000e+00 2.056671e-08 1.028336e-08
[86,] 1.000000e+00 1.267066e-08 6.335328e-09
[87,] 1.000000e+00 8.394382e-09 4.197191e-09
[88,] 1.000000e+00 5.998509e-09 2.999255e-09
[89,] 1.000000e+00 4.512210e-09 2.256105e-09
[90,] 1.000000e+00 3.343239e-09 1.671619e-09
[91,] 1.000000e+00 2.090059e-09 1.045029e-09
[92,] 1.000000e+00 7.176887e-10 3.588443e-10
[93,] 1.000000e+00 1.175831e-09 5.879153e-10
[94,] 1.000000e+00 1.719766e-09 8.598832e-10
[95,] 1.000000e+00 2.263904e-09 1.131952e-09
[96,] 1.000000e+00 2.793913e-09 1.396956e-09
[97,] 1.000000e+00 3.446347e-09 1.723174e-09
[98,] 1.000000e+00 4.552738e-09 2.276369e-09
[99,] 1.000000e+00 6.786329e-09 3.393164e-09
[100,] 1.000000e+00 1.160333e-08 5.801665e-09
[101,] 1.000000e+00 2.191558e-08 1.095779e-08
[102,] 1.000000e+00 3.982897e-08 1.991449e-08
[103,] 1.000000e+00 4.195052e-08 2.097526e-08
[104,] 1.000000e+00 3.612277e-112 1.806139e-112
[105,] 1.000000e+00 1.473645e-93 7.368224e-94
[106,] 1.000000e+00 2.779820e-79 1.389910e-79
[107,] 1.000000e+00 1.843712e-75 9.218560e-76
[108,] 1.000000e+00 1.689870e-52 8.449348e-53
[109,] 1.000000e+00 3.201643e-40 1.600822e-40
> postscript(file="/var/wessaorg/rcomp/tmp/10mpd1322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/29v0m1322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3b6qw1322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4o9u81322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5culw1322154290.ps",horizontal=F,onefile=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 = 118
Frequency = 1
1 2 3 4 5 6 7
-4.6034817 -4.5702879 -4.5370941 -4.5039003 -4.4707065 -4.4375127 -4.4043189
8 9 10 11 12 13 14
-4.3711251 -4.3379313 -4.3047375 -4.2715437 -4.2383499 -3.6034817 -3.5702879
15 16 17 18 19 20 21
-3.5370941 -3.5039003 -3.4707065 -3.4375127 -3.4043189 -3.3711251 -3.3379313
22 23 24 25 26 27 28
-3.3047375 -3.2715437 -3.2383499 -2.6034817 -2.5702879 -2.5370941 -2.5039003
29 30 31 32 33 34 35
-2.4707065 -2.4375127 -2.4043189 -2.3711251 -2.3379313 -2.3047375 -2.2715437
36 37 38 39 40 41 42
-2.2383499 -1.6034817 -1.5702879 -1.5370941 -1.5039003 -1.4707065 -1.4375127
43 44 45 46 47 48 49
-1.4043189 -1.3711251 -1.3379313 -1.3047375 -1.2715437 -1.2383499 -0.6034817
50 51 52 53 54 55 56
-0.5702879 -0.5370941 -0.5039003 -0.4707065 -0.4375127 -0.4043189 -0.3711251
57 58 59 60 61 62 63
-0.3379313 -0.3047375 -0.2715437 -0.2383499 0.3965183 0.4297121 0.4629059
64 65 66 67 68 69 70
0.4960997 0.5292935 0.5624873 0.5956811 0.6288749 0.6620687 0.6952625
71 72 73 74 75 76 77
0.7284563 0.7616501 1.3965183 1.4297121 1.4629059 1.4960997 1.5292935
78 79 80 81 82 83 84
1.5624873 1.5956811 1.6288749 1.6620687 1.6952625 1.7284563 1.7616501
85 86 87 88 89 90 91
2.3965183 2.4297121 2.4629059 2.4960997 2.5292935 2.5624873 2.5956811
92 93 94 95 96 97 98
2.6288749 2.6620687 2.6952625 2.7284563 2.7616501 3.3965183 3.4297121
99 100 101 102 103 104 105
3.4629059 3.4960997 3.5292935 3.5624873 3.5956811 3.6288749 3.6620687
106 107 108 109 110 111 112
3.6952625 3.7284563 3.7616501 4.3965183 4.4297121 4.4629059 4.4960997
113 114 115 116 117 118
4.5292935 4.5624873 4.5956811 4.6288749 4.6620687 4.6952625
> postscript(file="/var/wessaorg/rcomp/tmp/6jvvf1322154290.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 118
Frequency = 1
lag(myerror, k = 1) myerror
0 -4.6034817 NA
1 -4.5702879 -4.6034817
2 -4.5370941 -4.5702879
3 -4.5039003 -4.5370941
4 -4.4707065 -4.5039003
5 -4.4375127 -4.4707065
6 -4.4043189 -4.4375127
7 -4.3711251 -4.4043189
8 -4.3379313 -4.3711251
9 -4.3047375 -4.3379313
10 -4.2715437 -4.3047375
11 -4.2383499 -4.2715437
12 -3.6034817 -4.2383499
13 -3.5702879 -3.6034817
14 -3.5370941 -3.5702879
15 -3.5039003 -3.5370941
16 -3.4707065 -3.5039003
17 -3.4375127 -3.4707065
18 -3.4043189 -3.4375127
19 -3.3711251 -3.4043189
20 -3.3379313 -3.3711251
21 -3.3047375 -3.3379313
22 -3.2715437 -3.3047375
23 -3.2383499 -3.2715437
24 -2.6034817 -3.2383499
25 -2.5702879 -2.6034817
26 -2.5370941 -2.5702879
27 -2.5039003 -2.5370941
28 -2.4707065 -2.5039003
29 -2.4375127 -2.4707065
30 -2.4043189 -2.4375127
31 -2.3711251 -2.4043189
32 -2.3379313 -2.3711251
33 -2.3047375 -2.3379313
34 -2.2715437 -2.3047375
35 -2.2383499 -2.2715437
36 -1.6034817 -2.2383499
37 -1.5702879 -1.6034817
38 -1.5370941 -1.5702879
39 -1.5039003 -1.5370941
40 -1.4707065 -1.5039003
41 -1.4375127 -1.4707065
42 -1.4043189 -1.4375127
43 -1.3711251 -1.4043189
44 -1.3379313 -1.3711251
45 -1.3047375 -1.3379313
46 -1.2715437 -1.3047375
47 -1.2383499 -1.2715437
48 -0.6034817 -1.2383499
49 -0.5702879 -0.6034817
50 -0.5370941 -0.5702879
51 -0.5039003 -0.5370941
52 -0.4707065 -0.5039003
53 -0.4375127 -0.4707065
54 -0.4043189 -0.4375127
55 -0.3711251 -0.4043189
56 -0.3379313 -0.3711251
57 -0.3047375 -0.3379313
58 -0.2715437 -0.3047375
59 -0.2383499 -0.2715437
60 0.3965183 -0.2383499
61 0.4297121 0.3965183
62 0.4629059 0.4297121
63 0.4960997 0.4629059
64 0.5292935 0.4960997
65 0.5624873 0.5292935
66 0.5956811 0.5624873
67 0.6288749 0.5956811
68 0.6620687 0.6288749
69 0.6952625 0.6620687
70 0.7284563 0.6952625
71 0.7616501 0.7284563
72 1.3965183 0.7616501
73 1.4297121 1.3965183
74 1.4629059 1.4297121
75 1.4960997 1.4629059
76 1.5292935 1.4960997
77 1.5624873 1.5292935
78 1.5956811 1.5624873
79 1.6288749 1.5956811
80 1.6620687 1.6288749
81 1.6952625 1.6620687
82 1.7284563 1.6952625
83 1.7616501 1.7284563
84 2.3965183 1.7616501
85 2.4297121 2.3965183
86 2.4629059 2.4297121
87 2.4960997 2.4629059
88 2.5292935 2.4960997
89 2.5624873 2.5292935
90 2.5956811 2.5624873
91 2.6288749 2.5956811
92 2.6620687 2.6288749
93 2.6952625 2.6620687
94 2.7284563 2.6952625
95 2.7616501 2.7284563
96 3.3965183 2.7616501
97 3.4297121 3.3965183
98 3.4629059 3.4297121
99 3.4960997 3.4629059
100 3.5292935 3.4960997
101 3.5624873 3.5292935
102 3.5956811 3.5624873
103 3.6288749 3.5956811
104 3.6620687 3.6288749
105 3.6952625 3.6620687
106 3.7284563 3.6952625
107 3.7616501 3.7284563
108 4.3965183 3.7616501
109 4.4297121 4.3965183
110 4.4629059 4.4297121
111 4.4960997 4.4629059
112 4.5292935 4.4960997
113 4.5624873 4.5292935
114 4.5956811 4.5624873
115 4.6288749 4.5956811
116 4.6620687 4.6288749
117 4.6952625 4.6620687
118 NA 4.6952625
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4.5702879 -4.6034817
[2,] -4.5370941 -4.5702879
[3,] -4.5039003 -4.5370941
[4,] -4.4707065 -4.5039003
[5,] -4.4375127 -4.4707065
[6,] -4.4043189 -4.4375127
[7,] -4.3711251 -4.4043189
[8,] -4.3379313 -4.3711251
[9,] -4.3047375 -4.3379313
[10,] -4.2715437 -4.3047375
[11,] -4.2383499 -4.2715437
[12,] -3.6034817 -4.2383499
[13,] -3.5702879 -3.6034817
[14,] -3.5370941 -3.5702879
[15,] -3.5039003 -3.5370941
[16,] -3.4707065 -3.5039003
[17,] -3.4375127 -3.4707065
[18,] -3.4043189 -3.4375127
[19,] -3.3711251 -3.4043189
[20,] -3.3379313 -3.3711251
[21,] -3.3047375 -3.3379313
[22,] -3.2715437 -3.3047375
[23,] -3.2383499 -3.2715437
[24,] -2.6034817 -3.2383499
[25,] -2.5702879 -2.6034817
[26,] -2.5370941 -2.5702879
[27,] -2.5039003 -2.5370941
[28,] -2.4707065 -2.5039003
[29,] -2.4375127 -2.4707065
[30,] -2.4043189 -2.4375127
[31,] -2.3711251 -2.4043189
[32,] -2.3379313 -2.3711251
[33,] -2.3047375 -2.3379313
[34,] -2.2715437 -2.3047375
[35,] -2.2383499 -2.2715437
[36,] -1.6034817 -2.2383499
[37,] -1.5702879 -1.6034817
[38,] -1.5370941 -1.5702879
[39,] -1.5039003 -1.5370941
[40,] -1.4707065 -1.5039003
[41,] -1.4375127 -1.4707065
[42,] -1.4043189 -1.4375127
[43,] -1.3711251 -1.4043189
[44,] -1.3379313 -1.3711251
[45,] -1.3047375 -1.3379313
[46,] -1.2715437 -1.3047375
[47,] -1.2383499 -1.2715437
[48,] -0.6034817 -1.2383499
[49,] -0.5702879 -0.6034817
[50,] -0.5370941 -0.5702879
[51,] -0.5039003 -0.5370941
[52,] -0.4707065 -0.5039003
[53,] -0.4375127 -0.4707065
[54,] -0.4043189 -0.4375127
[55,] -0.3711251 -0.4043189
[56,] -0.3379313 -0.3711251
[57,] -0.3047375 -0.3379313
[58,] -0.2715437 -0.3047375
[59,] -0.2383499 -0.2715437
[60,] 0.3965183 -0.2383499
[61,] 0.4297121 0.3965183
[62,] 0.4629059 0.4297121
[63,] 0.4960997 0.4629059
[64,] 0.5292935 0.4960997
[65,] 0.5624873 0.5292935
[66,] 0.5956811 0.5624873
[67,] 0.6288749 0.5956811
[68,] 0.6620687 0.6288749
[69,] 0.6952625 0.6620687
[70,] 0.7284563 0.6952625
[71,] 0.7616501 0.7284563
[72,] 1.3965183 0.7616501
[73,] 1.4297121 1.3965183
[74,] 1.4629059 1.4297121
[75,] 1.4960997 1.4629059
[76,] 1.5292935 1.4960997
[77,] 1.5624873 1.5292935
[78,] 1.5956811 1.5624873
[79,] 1.6288749 1.5956811
[80,] 1.6620687 1.6288749
[81,] 1.6952625 1.6620687
[82,] 1.7284563 1.6952625
[83,] 1.7616501 1.7284563
[84,] 2.3965183 1.7616501
[85,] 2.4297121 2.3965183
[86,] 2.4629059 2.4297121
[87,] 2.4960997 2.4629059
[88,] 2.5292935 2.4960997
[89,] 2.5624873 2.5292935
[90,] 2.5956811 2.5624873
[91,] 2.6288749 2.5956811
[92,] 2.6620687 2.6288749
[93,] 2.6952625 2.6620687
[94,] 2.7284563 2.6952625
[95,] 2.7616501 2.7284563
[96,] 3.3965183 2.7616501
[97,] 3.4297121 3.3965183
[98,] 3.4629059 3.4297121
[99,] 3.4960997 3.4629059
[100,] 3.5292935 3.4960997
[101,] 3.5624873 3.5292935
[102,] 3.5956811 3.5624873
[103,] 3.6288749 3.5956811
[104,] 3.6620687 3.6288749
[105,] 3.6952625 3.6620687
[106,] 3.7284563 3.6952625
[107,] 3.7616501 3.7284563
[108,] 4.3965183 3.7616501
[109,] 4.4297121 4.3965183
[110,] 4.4629059 4.4297121
[111,] 4.4960997 4.4629059
[112,] 4.5292935 4.4960997
[113,] 4.5624873 4.5292935
[114,] 4.5956811 4.5624873
[115,] 4.6288749 4.5956811
[116,] 4.6620687 4.6288749
[117,] 4.6952625 4.6620687
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4.5702879 -4.6034817
2 -4.5370941 -4.5702879
3 -4.5039003 -4.5370941
4 -4.4707065 -4.5039003
5 -4.4375127 -4.4707065
6 -4.4043189 -4.4375127
7 -4.3711251 -4.4043189
8 -4.3379313 -4.3711251
9 -4.3047375 -4.3379313
10 -4.2715437 -4.3047375
11 -4.2383499 -4.2715437
12 -3.6034817 -4.2383499
13 -3.5702879 -3.6034817
14 -3.5370941 -3.5702879
15 -3.5039003 -3.5370941
16 -3.4707065 -3.5039003
17 -3.4375127 -3.4707065
18 -3.4043189 -3.4375127
19 -3.3711251 -3.4043189
20 -3.3379313 -3.3711251
21 -3.3047375 -3.3379313
22 -3.2715437 -3.3047375
23 -3.2383499 -3.2715437
24 -2.6034817 -3.2383499
25 -2.5702879 -2.6034817
26 -2.5370941 -2.5702879
27 -2.5039003 -2.5370941
28 -2.4707065 -2.5039003
29 -2.4375127 -2.4707065
30 -2.4043189 -2.4375127
31 -2.3711251 -2.4043189
32 -2.3379313 -2.3711251
33 -2.3047375 -2.3379313
34 -2.2715437 -2.3047375
35 -2.2383499 -2.2715437
36 -1.6034817 -2.2383499
37 -1.5702879 -1.6034817
38 -1.5370941 -1.5702879
39 -1.5039003 -1.5370941
40 -1.4707065 -1.5039003
41 -1.4375127 -1.4707065
42 -1.4043189 -1.4375127
43 -1.3711251 -1.4043189
44 -1.3379313 -1.3711251
45 -1.3047375 -1.3379313
46 -1.2715437 -1.3047375
47 -1.2383499 -1.2715437
48 -0.6034817 -1.2383499
49 -0.5702879 -0.6034817
50 -0.5370941 -0.5702879
51 -0.5039003 -0.5370941
52 -0.4707065 -0.5039003
53 -0.4375127 -0.4707065
54 -0.4043189 -0.4375127
55 -0.3711251 -0.4043189
56 -0.3379313 -0.3711251
57 -0.3047375 -0.3379313
58 -0.2715437 -0.3047375
59 -0.2383499 -0.2715437
60 0.3965183 -0.2383499
61 0.4297121 0.3965183
62 0.4629059 0.4297121
63 0.4960997 0.4629059
64 0.5292935 0.4960997
65 0.5624873 0.5292935
66 0.5956811 0.5624873
67 0.6288749 0.5956811
68 0.6620687 0.6288749
69 0.6952625 0.6620687
70 0.7284563 0.6952625
71 0.7616501 0.7284563
72 1.3965183 0.7616501
73 1.4297121 1.3965183
74 1.4629059 1.4297121
75 1.4960997 1.4629059
76 1.5292935 1.4960997
77 1.5624873 1.5292935
78 1.5956811 1.5624873
79 1.6288749 1.5956811
80 1.6620687 1.6288749
81 1.6952625 1.6620687
82 1.7284563 1.6952625
83 1.7616501 1.7284563
84 2.3965183 1.7616501
85 2.4297121 2.3965183
86 2.4629059 2.4297121
87 2.4960997 2.4629059
88 2.5292935 2.4960997
89 2.5624873 2.5292935
90 2.5956811 2.5624873
91 2.6288749 2.5956811
92 2.6620687 2.6288749
93 2.6952625 2.6620687
94 2.7284563 2.6952625
95 2.7616501 2.7284563
96 3.3965183 2.7616501
97 3.4297121 3.3965183
98 3.4629059 3.4297121
99 3.4960997 3.4629059
100 3.5292935 3.4960997
101 3.5624873 3.5292935
102 3.5956811 3.5624873
103 3.6288749 3.5956811
104 3.6620687 3.6288749
105 3.6952625 3.6620687
106 3.7284563 3.6952625
107 3.7616501 3.7284563
108 4.3965183 3.7616501
109 4.4297121 4.3965183
110 4.4629059 4.4297121
111 4.4960997 4.4629059
112 4.5292935 4.4960997
113 4.5624873 4.5292935
114 4.5956811 4.5624873
115 4.6288749 4.5956811
116 4.6620687 4.6288749
117 4.6952625 4.6620687
> 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/wessaorg/rcomp/tmp/7lzrk1322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/83y441322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9luhu1322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10qhh11322154290.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11o8yc1322154290.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/wessaorg/rcomp/tmp/12islu1322154290.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/wessaorg/rcomp/tmp/13n7r61322154290.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/wessaorg/rcomp/tmp/143nlo1322154290.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/wessaorg/rcomp/tmp/15md8q1322154290.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/wessaorg/rcomp/tmp/16cdux1322154290.tab")
+ }
>
> try(system("convert tmp/10mpd1322154290.ps tmp/10mpd1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/29v0m1322154290.ps tmp/29v0m1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/3b6qw1322154290.ps tmp/3b6qw1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/4o9u81322154290.ps tmp/4o9u81322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/5culw1322154290.ps tmp/5culw1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jvvf1322154290.ps tmp/6jvvf1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lzrk1322154290.ps tmp/7lzrk1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/83y441322154290.ps tmp/83y441322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/9luhu1322154290.ps tmp/9luhu1322154290.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qhh11322154290.ps tmp/10qhh11322154290.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.930 0.510 4.669