R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(7116,6927,6731,6850,6766,6979,7149,7067,7170,7237,7240,7645,7678,7491,7816,7631,8395,8578,8950,9450,9501,10083,10544,11299,12049,12860,13389,13796,14505,14727,14646,14861,15012,15421,15227,15124,14953,15039,15128,15221,14876,14517,14609,14735,14574,14636,15104,14393,13919,13751,13628,13792,13892,14024,13908,13920,13897,13759,13323,13097,12758,12806,12673,12500,12720,12749,12794,12544,12088,12258),dim=c(1,70),dimnames=list(c('Nojob'),1:70))
> y <- array(NA,dim=c(1,70),dimnames=list(c('Nojob'),1:70))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> par3 <- 'Linear Trend'
> par2 <- 'Include Monthly Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Nojob M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 7116 1 0 0 0 0 0 0 0 0 0 0 1
2 6927 0 1 0 0 0 0 0 0 0 0 0 2
3 6731 0 0 1 0 0 0 0 0 0 0 0 3
4 6850 0 0 0 1 0 0 0 0 0 0 0 4
5 6766 0 0 0 0 1 0 0 0 0 0 0 5
6 6979 0 0 0 0 0 1 0 0 0 0 0 6
7 7149 0 0 0 0 0 0 1 0 0 0 0 7
8 7067 0 0 0 0 0 0 0 1 0 0 0 8
9 7170 0 0 0 0 0 0 0 0 1 0 0 9
10 7237 0 0 0 0 0 0 0 0 0 1 0 10
11 7240 0 0 0 0 0 0 0 0 0 0 1 11
12 7645 0 0 0 0 0 0 0 0 0 0 0 12
13 7678 1 0 0 0 0 0 0 0 0 0 0 13
14 7491 0 1 0 0 0 0 0 0 0 0 0 14
15 7816 0 0 1 0 0 0 0 0 0 0 0 15
16 7631 0 0 0 1 0 0 0 0 0 0 0 16
17 8395 0 0 0 0 1 0 0 0 0 0 0 17
18 8578 0 0 0 0 0 1 0 0 0 0 0 18
19 8950 0 0 0 0 0 0 1 0 0 0 0 19
20 9450 0 0 0 0 0 0 0 1 0 0 0 20
21 9501 0 0 0 0 0 0 0 0 1 0 0 21
22 10083 0 0 0 0 0 0 0 0 0 1 0 22
23 10544 0 0 0 0 0 0 0 0 0 0 1 23
24 11299 0 0 0 0 0 0 0 0 0 0 0 24
25 12049 1 0 0 0 0 0 0 0 0 0 0 25
26 12860 0 1 0 0 0 0 0 0 0 0 0 26
27 13389 0 0 1 0 0 0 0 0 0 0 0 27
28 13796 0 0 0 1 0 0 0 0 0 0 0 28
29 14505 0 0 0 0 1 0 0 0 0 0 0 29
30 14727 0 0 0 0 0 1 0 0 0 0 0 30
31 14646 0 0 0 0 0 0 1 0 0 0 0 31
32 14861 0 0 0 0 0 0 0 1 0 0 0 32
33 15012 0 0 0 0 0 0 0 0 1 0 0 33
34 15421 0 0 0 0 0 0 0 0 0 1 0 34
35 15227 0 0 0 0 0 0 0 0 0 0 1 35
36 15124 0 0 0 0 0 0 0 0 0 0 0 36
37 14953 1 0 0 0 0 0 0 0 0 0 0 37
38 15039 0 1 0 0 0 0 0 0 0 0 0 38
39 15128 0 0 1 0 0 0 0 0 0 0 0 39
40 15221 0 0 0 1 0 0 0 0 0 0 0 40
41 14876 0 0 0 0 1 0 0 0 0 0 0 41
42 14517 0 0 0 0 0 1 0 0 0 0 0 42
43 14609 0 0 0 0 0 0 1 0 0 0 0 43
44 14735 0 0 0 0 0 0 0 1 0 0 0 44
45 14574 0 0 0 0 0 0 0 0 1 0 0 45
46 14636 0 0 0 0 0 0 0 0 0 1 0 46
47 15104 0 0 0 0 0 0 0 0 0 0 1 47
48 14393 0 0 0 0 0 0 0 0 0 0 0 48
49 13919 1 0 0 0 0 0 0 0 0 0 0 49
50 13751 0 1 0 0 0 0 0 0 0 0 0 50
51 13628 0 0 1 0 0 0 0 0 0 0 0 51
52 13792 0 0 0 1 0 0 0 0 0 0 0 52
53 13892 0 0 0 0 1 0 0 0 0 0 0 53
54 14024 0 0 0 0 0 1 0 0 0 0 0 54
55 13908 0 0 0 0 0 0 1 0 0 0 0 55
56 13920 0 0 0 0 0 0 0 1 0 0 0 56
57 13897 0 0 0 0 0 0 0 0 1 0 0 57
58 13759 0 0 0 0 0 0 0 0 0 1 0 58
59 13323 0 0 0 0 0 0 0 0 0 0 1 59
60 13097 0 0 0 0 0 0 0 0 0 0 0 60
61 12758 1 0 0 0 0 0 0 0 0 0 0 61
62 12806 0 1 0 0 0 0 0 0 0 0 0 62
63 12673 0 0 1 0 0 0 0 0 0 0 0 63
64 12500 0 0 0 1 0 0 0 0 0 0 0 64
65 12720 0 0 0 0 1 0 0 0 0 0 0 65
66 12749 0 0 0 0 0 1 0 0 0 0 0 66
67 12794 0 0 0 0 0 0 1 0 0 0 0 67
68 12544 0 0 0 0 0 0 0 1 0 0 0 68
69 12088 0 0 0 0 0 0 0 0 1 0 0 69
70 12258 0 0 0 0 0 0 0 0 0 1 0 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
8413.82 -358.08 -399.51 -425.95 -463.39 -344.33
M6 M7 M8 M9 M10 M11
-382.60 -410.54 -431.98 -596.08 -512.35 84.27
t
108.27
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3222.5 -1784.2 -311.3 1947.1 3838.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8413.82 1133.09 7.426 6.21e-10 ***
M1 -358.08 1384.59 -0.259 0.797
M2 -399.51 1383.99 -0.289 0.774
M3 -425.95 1383.52 -0.308 0.759
M4 -463.39 1383.18 -0.335 0.739
M5 -344.33 1382.98 -0.249 0.804
M6 -382.60 1382.92 -0.277 0.783
M7 -410.54 1382.98 -0.297 0.768
M8 -431.98 1383.18 -0.312 0.756
M9 -596.08 1383.52 -0.431 0.668
M10 -512.35 1383.99 -0.370 0.713
M11 84.27 1444.47 0.058 0.954
t 108.27 13.63 7.944 8.51e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2284 on 57 degrees of freedom
Multiple R-squared: 0.5302, Adjusted R-squared: 0.4313
F-statistic: 5.36 on 12 and 57 DF, p-value: 5.63e-06
> 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.0009860866 1.972173e-03 9.990139e-01
[2,] 0.0016623066 3.324613e-03 9.983377e-01
[3,] 0.0008547023 1.709405e-03 9.991453e-01
[4,] 0.0006309843 1.261969e-03 9.993690e-01
[5,] 0.0015401412 3.080282e-03 9.984599e-01
[6,] 0.0025513202 5.102640e-03 9.974487e-01
[7,] 0.0104879759 2.097595e-02 9.895120e-01
[8,] 0.0802429928 1.604860e-01 9.197570e-01
[9,] 0.3848998090 7.697996e-01 6.151002e-01
[10,] 0.7938753968 4.122492e-01 2.061246e-01
[11,] 0.9862801117 2.743978e-02 1.371989e-02
[12,] 0.9996802335 6.395331e-04 3.197665e-04
[13,] 0.9999973157 5.368634e-06 2.684317e-06
[14,] 0.9999998933 2.133418e-07 1.066709e-07
[15,] 0.9999999873 2.535909e-08 1.267955e-08
[16,] 0.9999999992 1.551141e-09 7.755704e-10
[17,] 0.9999999999 1.380231e-10 6.901155e-11
[18,] 1.0000000000 5.670396e-11 2.835198e-11
[19,] 0.9999999999 1.090208e-10 5.451038e-11
[20,] 1.0000000000 5.351083e-11 2.675541e-11
[21,] 0.9999999999 1.617885e-10 8.089424e-11
[22,] 0.9999999997 6.156680e-10 3.078340e-10
[23,] 0.9999999988 2.496934e-09 1.248467e-09
[24,] 0.9999999966 6.838701e-09 3.419350e-09
[25,] 0.9999999939 1.228199e-08 6.140996e-09
[26,] 0.9999999809 3.814056e-08 1.907028e-08
[27,] 0.9999999906 1.870346e-08 9.351730e-09
[28,] 0.9999999938 1.242629e-08 6.213147e-09
[29,] 0.9999999875 2.491831e-08 1.245916e-08
[30,] 0.9999999711 5.786589e-08 2.893294e-08
[31,] 0.9999999356 1.287917e-07 6.439586e-08
[32,] 0.9999999319 1.361367e-07 6.806837e-08
[33,] 0.9999996490 7.019372e-07 3.509686e-07
[34,] 0.9999985748 2.850386e-06 1.425193e-06
[35,] 0.9999967093 6.581369e-06 3.290685e-06
[36,] 0.9999932393 1.352131e-05 6.760653e-06
[37,] 0.9999479953 1.040094e-04 5.200470e-05
[38,] 0.9996922215 6.155570e-04 3.077785e-04
[39,] 0.9976334477 4.733105e-03 2.366552e-03
> postscript(file="/var/fisher/rcomp/tmp/146wg1355847581.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/fisher/rcomp/tmp/2asxd1355847581.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/fisher/rcomp/tmp/329oi1355847581.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/fisher/rcomp/tmp/4399l1355847581.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/fisher/rcomp/tmp/518z21355847581.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 = 70
Frequency = 1
1 2 3 4 5 6
-1048.01923 -1303.85256 -1581.68590 -1533.51923 -1844.85256 -1701.85256
7 8 9 10 11 12
-1612.18590 -1781.01923 -1622.18590 -1747.18590 -2449.08205 -2068.08205
13 14 15 16 17 18
-1785.27821 -2039.11154 -1795.94487 -2051.77821 -1515.11154 -1402.11154
19 20 21 22 23 24
-1110.44487 -697.27821 -590.44487 -200.44487 -444.34103 286.65897
25 26 27 28 29 30
1286.46282 2030.62949 2477.79615 2813.96282 3295.62949 3447.62949
31 32 33 34 35 36
3286.29615 3414.46282 3621.29615 3838.29615 2939.40000 2812.40000
37 38 39 40 41 42
2891.20385 2910.37051 2917.53718 2939.70385 2367.37051 1938.37051
43 44 45 46 47 48
1950.03718 1989.20385 1884.03718 1754.03718 1517.14103 782.14103
49 50 51 52 53 54
557.94487 323.11154 118.27821 211.44487 84.11154 146.11154
55 56 57 58 59 60
-50.22179 -125.05513 -92.22179 -422.22179 -1563.11795 -1813.11795
61 62 63 64 65 66
-1902.31410 -1921.14744 -2135.98077 -2379.81410 -2387.14744 -2428.14744
67 68 69 70
-2463.48077 -2800.31410 -3200.48077 -3222.48077
> postscript(file="/var/fisher/rcomp/tmp/63fbi1355847581.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -1048.01923 NA
1 -1303.85256 -1048.01923
2 -1581.68590 -1303.85256
3 -1533.51923 -1581.68590
4 -1844.85256 -1533.51923
5 -1701.85256 -1844.85256
6 -1612.18590 -1701.85256
7 -1781.01923 -1612.18590
8 -1622.18590 -1781.01923
9 -1747.18590 -1622.18590
10 -2449.08205 -1747.18590
11 -2068.08205 -2449.08205
12 -1785.27821 -2068.08205
13 -2039.11154 -1785.27821
14 -1795.94487 -2039.11154
15 -2051.77821 -1795.94487
16 -1515.11154 -2051.77821
17 -1402.11154 -1515.11154
18 -1110.44487 -1402.11154
19 -697.27821 -1110.44487
20 -590.44487 -697.27821
21 -200.44487 -590.44487
22 -444.34103 -200.44487
23 286.65897 -444.34103
24 1286.46282 286.65897
25 2030.62949 1286.46282
26 2477.79615 2030.62949
27 2813.96282 2477.79615
28 3295.62949 2813.96282
29 3447.62949 3295.62949
30 3286.29615 3447.62949
31 3414.46282 3286.29615
32 3621.29615 3414.46282
33 3838.29615 3621.29615
34 2939.40000 3838.29615
35 2812.40000 2939.40000
36 2891.20385 2812.40000
37 2910.37051 2891.20385
38 2917.53718 2910.37051
39 2939.70385 2917.53718
40 2367.37051 2939.70385
41 1938.37051 2367.37051
42 1950.03718 1938.37051
43 1989.20385 1950.03718
44 1884.03718 1989.20385
45 1754.03718 1884.03718
46 1517.14103 1754.03718
47 782.14103 1517.14103
48 557.94487 782.14103
49 323.11154 557.94487
50 118.27821 323.11154
51 211.44487 118.27821
52 84.11154 211.44487
53 146.11154 84.11154
54 -50.22179 146.11154
55 -125.05513 -50.22179
56 -92.22179 -125.05513
57 -422.22179 -92.22179
58 -1563.11795 -422.22179
59 -1813.11795 -1563.11795
60 -1902.31410 -1813.11795
61 -1921.14744 -1902.31410
62 -2135.98077 -1921.14744
63 -2379.81410 -2135.98077
64 -2387.14744 -2379.81410
65 -2428.14744 -2387.14744
66 -2463.48077 -2428.14744
67 -2800.31410 -2463.48077
68 -3200.48077 -2800.31410
69 -3222.48077 -3200.48077
70 NA -3222.48077
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1303.85256 -1048.01923
[2,] -1581.68590 -1303.85256
[3,] -1533.51923 -1581.68590
[4,] -1844.85256 -1533.51923
[5,] -1701.85256 -1844.85256
[6,] -1612.18590 -1701.85256
[7,] -1781.01923 -1612.18590
[8,] -1622.18590 -1781.01923
[9,] -1747.18590 -1622.18590
[10,] -2449.08205 -1747.18590
[11,] -2068.08205 -2449.08205
[12,] -1785.27821 -2068.08205
[13,] -2039.11154 -1785.27821
[14,] -1795.94487 -2039.11154
[15,] -2051.77821 -1795.94487
[16,] -1515.11154 -2051.77821
[17,] -1402.11154 -1515.11154
[18,] -1110.44487 -1402.11154
[19,] -697.27821 -1110.44487
[20,] -590.44487 -697.27821
[21,] -200.44487 -590.44487
[22,] -444.34103 -200.44487
[23,] 286.65897 -444.34103
[24,] 1286.46282 286.65897
[25,] 2030.62949 1286.46282
[26,] 2477.79615 2030.62949
[27,] 2813.96282 2477.79615
[28,] 3295.62949 2813.96282
[29,] 3447.62949 3295.62949
[30,] 3286.29615 3447.62949
[31,] 3414.46282 3286.29615
[32,] 3621.29615 3414.46282
[33,] 3838.29615 3621.29615
[34,] 2939.40000 3838.29615
[35,] 2812.40000 2939.40000
[36,] 2891.20385 2812.40000
[37,] 2910.37051 2891.20385
[38,] 2917.53718 2910.37051
[39,] 2939.70385 2917.53718
[40,] 2367.37051 2939.70385
[41,] 1938.37051 2367.37051
[42,] 1950.03718 1938.37051
[43,] 1989.20385 1950.03718
[44,] 1884.03718 1989.20385
[45,] 1754.03718 1884.03718
[46,] 1517.14103 1754.03718
[47,] 782.14103 1517.14103
[48,] 557.94487 782.14103
[49,] 323.11154 557.94487
[50,] 118.27821 323.11154
[51,] 211.44487 118.27821
[52,] 84.11154 211.44487
[53,] 146.11154 84.11154
[54,] -50.22179 146.11154
[55,] -125.05513 -50.22179
[56,] -92.22179 -125.05513
[57,] -422.22179 -92.22179
[58,] -1563.11795 -422.22179
[59,] -1813.11795 -1563.11795
[60,] -1902.31410 -1813.11795
[61,] -1921.14744 -1902.31410
[62,] -2135.98077 -1921.14744
[63,] -2379.81410 -2135.98077
[64,] -2387.14744 -2379.81410
[65,] -2428.14744 -2387.14744
[66,] -2463.48077 -2428.14744
[67,] -2800.31410 -2463.48077
[68,] -3200.48077 -2800.31410
[69,] -3222.48077 -3200.48077
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1303.85256 -1048.01923
2 -1581.68590 -1303.85256
3 -1533.51923 -1581.68590
4 -1844.85256 -1533.51923
5 -1701.85256 -1844.85256
6 -1612.18590 -1701.85256
7 -1781.01923 -1612.18590
8 -1622.18590 -1781.01923
9 -1747.18590 -1622.18590
10 -2449.08205 -1747.18590
11 -2068.08205 -2449.08205
12 -1785.27821 -2068.08205
13 -2039.11154 -1785.27821
14 -1795.94487 -2039.11154
15 -2051.77821 -1795.94487
16 -1515.11154 -2051.77821
17 -1402.11154 -1515.11154
18 -1110.44487 -1402.11154
19 -697.27821 -1110.44487
20 -590.44487 -697.27821
21 -200.44487 -590.44487
22 -444.34103 -200.44487
23 286.65897 -444.34103
24 1286.46282 286.65897
25 2030.62949 1286.46282
26 2477.79615 2030.62949
27 2813.96282 2477.79615
28 3295.62949 2813.96282
29 3447.62949 3295.62949
30 3286.29615 3447.62949
31 3414.46282 3286.29615
32 3621.29615 3414.46282
33 3838.29615 3621.29615
34 2939.40000 3838.29615
35 2812.40000 2939.40000
36 2891.20385 2812.40000
37 2910.37051 2891.20385
38 2917.53718 2910.37051
39 2939.70385 2917.53718
40 2367.37051 2939.70385
41 1938.37051 2367.37051
42 1950.03718 1938.37051
43 1989.20385 1950.03718
44 1884.03718 1989.20385
45 1754.03718 1884.03718
46 1517.14103 1754.03718
47 782.14103 1517.14103
48 557.94487 782.14103
49 323.11154 557.94487
50 118.27821 323.11154
51 211.44487 118.27821
52 84.11154 211.44487
53 146.11154 84.11154
54 -50.22179 146.11154
55 -125.05513 -50.22179
56 -92.22179 -125.05513
57 -422.22179 -92.22179
58 -1563.11795 -422.22179
59 -1813.11795 -1563.11795
60 -1902.31410 -1813.11795
61 -1921.14744 -1902.31410
62 -2135.98077 -1921.14744
63 -2379.81410 -2135.98077
64 -2387.14744 -2379.81410
65 -2428.14744 -2387.14744
66 -2463.48077 -2428.14744
67 -2800.31410 -2463.48077
68 -3200.48077 -2800.31410
69 -3222.48077 -3200.48077
> 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/fisher/rcomp/tmp/75an71355847581.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/fisher/rcomp/tmp/81qvk1355847581.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/fisher/rcomp/tmp/96ks41355847581.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/fisher/rcomp/tmp/10g9lo1355847581.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11ps931355847581.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/fisher/rcomp/tmp/12ppn11355847581.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/fisher/rcomp/tmp/13gr611355847581.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/fisher/rcomp/tmp/14zifz1355847581.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/fisher/rcomp/tmp/15q7lw1355847581.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/fisher/rcomp/tmp/168tbr1355847581.tab")
+ }
>
> try(system("convert tmp/146wg1355847581.ps tmp/146wg1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/2asxd1355847581.ps tmp/2asxd1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/329oi1355847581.ps tmp/329oi1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/4399l1355847581.ps tmp/4399l1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/518z21355847581.ps tmp/518z21355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/63fbi1355847581.ps tmp/63fbi1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/75an71355847581.ps tmp/75an71355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/81qvk1355847581.ps tmp/81qvk1355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/96ks41355847581.ps tmp/96ks41355847581.png",intern=TRUE))
character(0)
> try(system("convert tmp/10g9lo1355847581.ps tmp/10g9lo1355847581.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.008 1.621 7.631