R version 3.6.2 (2019-12-12) -- "Dark and Stormy Night"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-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.
> source('/home/pw/wessanet/cretab')
>
>
> RC.capture <- function (expression, collapse = NULL) {
+ resultConn <- textConnection('RC.resultText', open = 'w', local=TRUE)
+ sink(resultConn)
+ on.exit(function() {
+ sink()
+ close(resultConn)
+ })
+ expression
+ on.exit(NULL)
+ sink()
+ close(resultConn)
+ return(paste(c(RC.resultText, ''), collapse = collapse, sep = ''))
+ }
> RC.texteval <- function (sourceText, collapse = NULL, echo = TRUE) {
+ sourceConn <- textConnection(sourceText, open = 'r')
+ on.exit(close(sourceConn))
+ result <- RC.capture(source(file = sourceConn, local = FALSE, echo = echo, print.eval = TRUE), collapse = collapse)
+ on.exit(NULL)
+ close(sourceConn)
+ res <- ''
+ for(i in 1:length(result)) {
+ if (result[i]!='') res <- paste(res,result[i],'
+ ',sep='')
+ }
+ return(res)
+ }
>
>
> myrfcuid = ''
>
> x <- array(list(.100,1.63,.070,0.91,.180,1.24,.030,1.39,-.077,0.50,.100,0.75,.123,0.23,.047,0.19,.063,0.40,-.043,0.15,-.070,1.25,-.083,1.42,.036,1.51,-.036,0.72,.020,0.59,.080,0.32,.003,0.54,.177,0.22,.280,0.06,.233,0.61,.227,0.31,.236,0.03,.250,-0.01,.137,-0.63,.337,-0.20,.193,1.47,-.017,1.46,.030,1.78,.050,1.86,.050,1.20,.150,1.00,.177,-1.26,.027,-0.37,.176,-0.30,.017,1.33,.190,-0.10,.060,0.70,-.010,1.03,.117,0.84,.043,1.30,.117,0.93,.113,0.97,.087,-.13,.073,0.80,.097,1.53,.056,1.37,.187,1.53,.083,1.47,.127,1.00,.050,1.06,.070,2.54,.143,2.66,.097,1.20,.200,0.94,.093,1.86,.167,3.00,.297,2.90,.290,1.84,.280,-.54,.266,0.50,.084,1.70,.200,2.40,.333,3.87,.367,2.93,.333,2.40,.333,3.17),dim=c(2,66),dimnames=list(c('^CPI','^M'),1:66))
> y <- array(NA,dim=c(2,66),dimnames=list(c('^CPI','^M'),1:66))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par6 = '12'
> par5 = '0'
> par4 = '0'
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par6 <- '12'
> par5 <- '0'
> par4 <- '0'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.2.327 (Fri, 21 Jul 2017 20:18:06 +0200)
> #Author: root
> #To cite this work: Wessa P., (2017), Multiple Regression (v1.0.48) in Free Statistics Software (v$_version), Office for Research Development and Education, URL https://www.wessa.net/rwasp_multipleregression.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
> library(car)
Loading required package: carData
> library(MASS)
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> mywarning <- ''
> par6 <- as.numeric(par6)
> if(is.na(par6)) {
+ par6 <- 12
+ mywarning = 'Warning: you did not specify the seasonality. The seasonal period was set to s = 12.'
+ }
> par1 <- as.numeric(par1)
> if(is.na(par1)) {
+ par1 <- 1
+ mywarning = 'Warning: you did not specify the column number of the endogenous series! The first column was selected by default.'
+ }
> if (par4=='') par4 <- 0
> par4 <- as.numeric(par4)
> if (!is.numeric(par4)) par4 <- 0
> if (par5=='') par5 <- 0
> par5 <- as.numeric(par5)
> if (!is.numeric(par5)) par5 <- 0
> x <- na.omit(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'){
+ (n <- n -1)
+ x2 <- array(0, dim=c(n,k), dimnames=list(1:n, paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par3 == 'Seasonal Differences (s)'){
+ (n <- n - par6)
+ x2 <- array(0, dim=c(n,k), dimnames=list(1:n, paste('(1-Bs)',colnames(x),sep='')))
+ for (i in 1:n) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+par6,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par3 == 'First and Seasonal Differences (s)'){
+ (n <- n -1)
+ x2 <- array(0, dim=c(n,k), dimnames=list(1:n, paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ (n <- n - par6)
+ x2 <- array(0, dim=c(n,k), dimnames=list(1:n, paste('(1-Bs)',colnames(x),sep='')))
+ for (i in 1:n) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+par6,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if(par4 > 0) {
+ x2 <- array(0, dim=c(n-par4,par4), dimnames=list(1:(n-par4), paste(colnames(x)[par1],'(t-',1:par4,')',sep='')))
+ for (i in 1:(n-par4)) {
+ for (j in 1:par4) {
+ x2[i,j] <- x[i+par4-j,par1]
+ }
+ }
+ x <- cbind(x[(par4+1):n,], x2)
+ n <- n - par4
+ }
> if(par5 > 0) {
+ x2 <- array(0, dim=c(n-par5*par6,par5), dimnames=list(1:(n-par5*par6), paste(colnames(x)[par1],'(t-',1:par5,'s)',sep='')))
+ for (i in 1:(n-par5*par6)) {
+ for (j in 1:par5) {
+ x2[i,j] <- x[i+par5*par6-j*par6,par1]
+ }
+ }
+ x <- cbind(x[(par5*par6+1):n,], x2)
+ n <- n - par5*par6
+ }
> if (par2 == 'Include Seasonal Dummies'){
+ x2 <- array(0, dim=c(n,par6-1), dimnames=list(1:n, paste('M', seq(1:(par6-1)), sep ='')))
+ for (i in 1:(par6-1)){
+ x2[seq(i,n,par6),i] <- 1
+ }
+ x <- cbind(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[n,]))
[1] 2
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> print(x)
^CPI ^M
1 0.100 1.63
2 0.070 0.91
3 0.180 1.24
4 0.030 1.39
5 -0.077 0.50
6 0.100 0.75
7 0.123 0.23
8 0.047 0.19
9 0.063 0.40
10 -0.043 0.15
11 -0.070 1.25
12 -0.083 1.42
13 0.036 1.51
14 -0.036 0.72
15 0.020 0.59
16 0.080 0.32
17 0.003 0.54
18 0.177 0.22
19 0.280 0.06
20 0.233 0.61
21 0.227 0.31
22 0.236 0.03
23 0.250 -0.01
24 0.137 -0.63
25 0.337 -0.20
26 0.193 1.47
27 -0.017 1.46
28 0.030 1.78
29 0.050 1.86
30 0.050 1.20
31 0.150 1.00
32 0.177 -1.26
33 0.027 -0.37
34 0.176 -0.30
35 0.017 1.33
36 0.190 -0.10
37 0.060 0.70
38 -0.010 1.03
39 0.117 0.84
40 0.043 1.30
41 0.117 0.93
42 0.113 0.97
43 0.087 -0.13
44 0.073 0.80
45 0.097 1.53
46 0.056 1.37
47 0.187 1.53
48 0.083 1.47
49 0.127 1.00
50 0.050 1.06
51 0.070 2.54
52 0.143 2.66
53 0.097 1.20
54 0.200 0.94
55 0.093 1.86
56 0.167 3.00
57 0.297 2.90
58 0.290 1.84
59 0.280 -0.54
60 0.266 0.50
61 0.084 1.70
62 0.200 2.40
63 0.333 3.87
64 0.367 2.93
65 0.333 2.40
66 0.333 3.17
> (k <- length(x[n,]))
[1] 2
> head(x)
^CPI ^M
1 0.100 1.63
2 0.070 0.91
3 0.180 1.24
4 0.030 1.39
5 -0.077 0.50
6 0.100 0.75
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `^M`
0.10473 0.01863
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.21419 -0.07644 -0.01701 0.07753 0.23600
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.10473 0.01982 5.285 1.62e-06 ***
`^M` 0.01863 0.01366 1.365 0.177
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1099 on 64 degrees of freedom
Multiple R-squared: 0.02827, Adjusted R-squared: 0.01309
F-statistic: 1.862 on 1 and 64 DF, p-value: 0.1772
> 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.3535778 0.7071555 0.6464222
[2,] 0.2912378 0.5824756 0.7087622
[3,] 0.2977934 0.5955869 0.7022066
[4,] 0.1884911 0.3769821 0.8115089
[5,] 0.1122935 0.2245871 0.8877065
[6,] 0.1043492 0.2086984 0.8956508
[7,] 0.2058694 0.4117388 0.7941306
[8,] 0.3120451 0.6240901 0.6879549
[9,] 0.2412208 0.4824417 0.7587792
[10,] 0.2334743 0.4669487 0.7665257
[11,] 0.1821430 0.3642860 0.8178570
[12,] 0.1425040 0.2850079 0.8574960
[13,] 0.1169155 0.2338310 0.8830845
[14,] 0.1580647 0.3161295 0.8419353
[15,] 0.3677732 0.7355463 0.6322268
[16,] 0.4663677 0.9327355 0.5336323
[17,] 0.5023196 0.9953608 0.4976804
[18,] 0.5148541 0.9702918 0.4851459
[19,] 0.5347350 0.9305299 0.4652650
[20,] 0.4760997 0.9521994 0.5239003
[21,] 0.6724280 0.6551441 0.3275720
[22,] 0.7169289 0.5661422 0.2830711
[23,] 0.7276534 0.5446933 0.2723466
[24,] 0.7134171 0.5731658 0.2865829
[25,] 0.6967888 0.6064224 0.3032112
[26,] 0.6605808 0.6788383 0.3394192
[27,] 0.6136560 0.7726880 0.3863440
[28,] 0.6082178 0.7835644 0.3917822
[29,] 0.6079726 0.7840549 0.3920274
[30,] 0.5728460 0.8543080 0.4271540
[31,] 0.5729943 0.8540115 0.4270057
[32,] 0.5566577 0.8866847 0.4433423
[33,] 0.5006831 0.9986338 0.4993169
[34,] 0.5451144 0.9097711 0.4548856
[35,] 0.4771965 0.9543929 0.5228035
[36,] 0.4647299 0.9294599 0.5352701
[37,] 0.4003542 0.8007084 0.5996458
[38,] 0.3395701 0.6791401 0.6604299
[39,] 0.2826565 0.5653130 0.7173435
[40,] 0.2430592 0.4861184 0.7569408
[41,] 0.2178992 0.4357983 0.7821008
[42,] 0.2200782 0.4401564 0.7799218
[43,] 0.2046551 0.4093102 0.7953449
[44,] 0.1944529 0.3889058 0.8055471
[45,] 0.1570554 0.3141107 0.8429446
[46,] 0.1866247 0.3732493 0.8133753
[47,] 0.2755405 0.5510810 0.7244595
[48,] 0.3110149 0.6220298 0.6889851
[49,] 0.3523135 0.7046269 0.6476865
[50,] 0.2957934 0.5915868 0.7042066
[51,] 0.4493170 0.8986340 0.5506830
[52,] 0.5526876 0.8946247 0.4473124
[53,] 0.5260960 0.9478081 0.4739040
[54,] 0.4590946 0.9181892 0.5409054
[55,] 0.4594017 0.9188034 0.5405983
[56,] 0.6237259 0.7525481 0.3762741
[57,] 0.7120043 0.5759914 0.2879957
> postscript(file="/home/pw/wessanet/rcomp/tmp/1zr511582831289.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="/home/pw/wessanet/rcomp/tmp/21nk91582831289.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="/home/pw/wessanet/rcomp/tmp/3cmva1582831289.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> sresid <- studres(mylm)
> hist(sresid, freq=FALSE, main='Distribution of Studentized Residuals')
> xfit<-seq(min(sresid),max(sresid),length=40)
> yfit<-dnorm(xfit)
> lines(xfit, yfit)
> grid()
> dev.off()
null device
1
> postscript(file="/home/pw/wessanet/rcomp/tmp/4bcj21582831289.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="/home/pw/wessanet/rcomp/tmp/5ge8k1582831289.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqPlot(mylm, main='QQ Plot')
[1] 12 25
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 66
Frequency = 1
1 2 3 4 5 6
-0.035103287 -0.051686956 0.052163892 -0.100631176 -0.191047102 -0.018705550
7 8 9 10 11 12
0.013984022 -0.061270626 -0.049183723 -0.150525275 -0.198022446 -0.214190190
13 14 15 16 17 18
-0.096867232 -0.154146536 -0.095724143 -0.030693019 -0.111792453 0.068170360
19 20 21 22 23 24
0.174151767 0.116903181 0.116493319 0.130710780 0.145456132 0.044009083
25 26 27 28 29 30
0.235996553 0.060878120 -0.148935542 -0.107898355 -0.089389059 -0.077090756
31 32 33 34 35 36
0.026636002 0.095748372 -0.070835703 0.076859932 -0.112513149 0.087133173
37 38 39 40 41 42
-0.057773860 -0.133923011 -0.003382591 -0.085954135 -0.005059632 -0.009804984
43 44 45 46 47 48
-0.015307813 -0.046637239 -0.036239907 -0.074258501 0.053760093 -0.049121880
49 50 51 52 53 54
0.003636002 -0.074482025 -0.082060037 -0.011296092 -0.030090756 0.077754030
55 56 57 58 59 60
-0.046389059 0.006368419 0.138231798 0.150983617 0.185332042 0.151952898
61 62 63 64 65 66
-0.052407652 0.050548694 0.156157020 0.207672784 0.183548694 0.169200674
> postscript(file="/home/pw/wessanet/rcomp/tmp/6e2wo1582831289.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 = 66
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.035103287 NA
1 -0.051686956 -0.035103287
2 0.052163892 -0.051686956
3 -0.100631176 0.052163892
4 -0.191047102 -0.100631176
5 -0.018705550 -0.191047102
6 0.013984022 -0.018705550
7 -0.061270626 0.013984022
8 -0.049183723 -0.061270626
9 -0.150525275 -0.049183723
10 -0.198022446 -0.150525275
11 -0.214190190 -0.198022446
12 -0.096867232 -0.214190190
13 -0.154146536 -0.096867232
14 -0.095724143 -0.154146536
15 -0.030693019 -0.095724143
16 -0.111792453 -0.030693019
17 0.068170360 -0.111792453
18 0.174151767 0.068170360
19 0.116903181 0.174151767
20 0.116493319 0.116903181
21 0.130710780 0.116493319
22 0.145456132 0.130710780
23 0.044009083 0.145456132
24 0.235996553 0.044009083
25 0.060878120 0.235996553
26 -0.148935542 0.060878120
27 -0.107898355 -0.148935542
28 -0.089389059 -0.107898355
29 -0.077090756 -0.089389059
30 0.026636002 -0.077090756
31 0.095748372 0.026636002
32 -0.070835703 0.095748372
33 0.076859932 -0.070835703
34 -0.112513149 0.076859932
35 0.087133173 -0.112513149
36 -0.057773860 0.087133173
37 -0.133923011 -0.057773860
38 -0.003382591 -0.133923011
39 -0.085954135 -0.003382591
40 -0.005059632 -0.085954135
41 -0.009804984 -0.005059632
42 -0.015307813 -0.009804984
43 -0.046637239 -0.015307813
44 -0.036239907 -0.046637239
45 -0.074258501 -0.036239907
46 0.053760093 -0.074258501
47 -0.049121880 0.053760093
48 0.003636002 -0.049121880
49 -0.074482025 0.003636002
50 -0.082060037 -0.074482025
51 -0.011296092 -0.082060037
52 -0.030090756 -0.011296092
53 0.077754030 -0.030090756
54 -0.046389059 0.077754030
55 0.006368419 -0.046389059
56 0.138231798 0.006368419
57 0.150983617 0.138231798
58 0.185332042 0.150983617
59 0.151952898 0.185332042
60 -0.052407652 0.151952898
61 0.050548694 -0.052407652
62 0.156157020 0.050548694
63 0.207672784 0.156157020
64 0.183548694 0.207672784
65 0.169200674 0.183548694
66 NA 0.169200674
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.051686956 -0.035103287
[2,] 0.052163892 -0.051686956
[3,] -0.100631176 0.052163892
[4,] -0.191047102 -0.100631176
[5,] -0.018705550 -0.191047102
[6,] 0.013984022 -0.018705550
[7,] -0.061270626 0.013984022
[8,] -0.049183723 -0.061270626
[9,] -0.150525275 -0.049183723
[10,] -0.198022446 -0.150525275
[11,] -0.214190190 -0.198022446
[12,] -0.096867232 -0.214190190
[13,] -0.154146536 -0.096867232
[14,] -0.095724143 -0.154146536
[15,] -0.030693019 -0.095724143
[16,] -0.111792453 -0.030693019
[17,] 0.068170360 -0.111792453
[18,] 0.174151767 0.068170360
[19,] 0.116903181 0.174151767
[20,] 0.116493319 0.116903181
[21,] 0.130710780 0.116493319
[22,] 0.145456132 0.130710780
[23,] 0.044009083 0.145456132
[24,] 0.235996553 0.044009083
[25,] 0.060878120 0.235996553
[26,] -0.148935542 0.060878120
[27,] -0.107898355 -0.148935542
[28,] -0.089389059 -0.107898355
[29,] -0.077090756 -0.089389059
[30,] 0.026636002 -0.077090756
[31,] 0.095748372 0.026636002
[32,] -0.070835703 0.095748372
[33,] 0.076859932 -0.070835703
[34,] -0.112513149 0.076859932
[35,] 0.087133173 -0.112513149
[36,] -0.057773860 0.087133173
[37,] -0.133923011 -0.057773860
[38,] -0.003382591 -0.133923011
[39,] -0.085954135 -0.003382591
[40,] -0.005059632 -0.085954135
[41,] -0.009804984 -0.005059632
[42,] -0.015307813 -0.009804984
[43,] -0.046637239 -0.015307813
[44,] -0.036239907 -0.046637239
[45,] -0.074258501 -0.036239907
[46,] 0.053760093 -0.074258501
[47,] -0.049121880 0.053760093
[48,] 0.003636002 -0.049121880
[49,] -0.074482025 0.003636002
[50,] -0.082060037 -0.074482025
[51,] -0.011296092 -0.082060037
[52,] -0.030090756 -0.011296092
[53,] 0.077754030 -0.030090756
[54,] -0.046389059 0.077754030
[55,] 0.006368419 -0.046389059
[56,] 0.138231798 0.006368419
[57,] 0.150983617 0.138231798
[58,] 0.185332042 0.150983617
[59,] 0.151952898 0.185332042
[60,] -0.052407652 0.151952898
[61,] 0.050548694 -0.052407652
[62,] 0.156157020 0.050548694
[63,] 0.207672784 0.156157020
[64,] 0.183548694 0.207672784
[65,] 0.169200674 0.183548694
> z <- as.data.frame(dum1)
> print(z)
lag(myerror, k = 1) myerror
1 -0.051686956 -0.035103287
2 0.052163892 -0.051686956
3 -0.100631176 0.052163892
4 -0.191047102 -0.100631176
5 -0.018705550 -0.191047102
6 0.013984022 -0.018705550
7 -0.061270626 0.013984022
8 -0.049183723 -0.061270626
9 -0.150525275 -0.049183723
10 -0.198022446 -0.150525275
11 -0.214190190 -0.198022446
12 -0.096867232 -0.214190190
13 -0.154146536 -0.096867232
14 -0.095724143 -0.154146536
15 -0.030693019 -0.095724143
16 -0.111792453 -0.030693019
17 0.068170360 -0.111792453
18 0.174151767 0.068170360
19 0.116903181 0.174151767
20 0.116493319 0.116903181
21 0.130710780 0.116493319
22 0.145456132 0.130710780
23 0.044009083 0.145456132
24 0.235996553 0.044009083
25 0.060878120 0.235996553
26 -0.148935542 0.060878120
27 -0.107898355 -0.148935542
28 -0.089389059 -0.107898355
29 -0.077090756 -0.089389059
30 0.026636002 -0.077090756
31 0.095748372 0.026636002
32 -0.070835703 0.095748372
33 0.076859932 -0.070835703
34 -0.112513149 0.076859932
35 0.087133173 -0.112513149
36 -0.057773860 0.087133173
37 -0.133923011 -0.057773860
38 -0.003382591 -0.133923011
39 -0.085954135 -0.003382591
40 -0.005059632 -0.085954135
41 -0.009804984 -0.005059632
42 -0.015307813 -0.009804984
43 -0.046637239 -0.015307813
44 -0.036239907 -0.046637239
45 -0.074258501 -0.036239907
46 0.053760093 -0.074258501
47 -0.049121880 0.053760093
48 0.003636002 -0.049121880
49 -0.074482025 0.003636002
50 -0.082060037 -0.074482025
51 -0.011296092 -0.082060037
52 -0.030090756 -0.011296092
53 0.077754030 -0.030090756
54 -0.046389059 0.077754030
55 0.006368419 -0.046389059
56 0.138231798 0.006368419
57 0.150983617 0.138231798
58 0.185332042 0.150983617
59 0.151952898 0.185332042
60 -0.052407652 0.151952898
61 0.050548694 -0.052407652
62 0.156157020 0.050548694
63 0.207672784 0.156157020
64 0.183548694 0.207672784
65 0.169200674 0.183548694
> 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="/home/pw/wessanet/rcomp/tmp/7j3au1582831289.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="/home/pw/wessanet/rcomp/tmp/8kg0i1582831289.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="/home/pw/wessanet/rcomp/tmp/96oss1582831289.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="/home/pw/wessanet/rcomp/tmp/106e0m1582831289.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
>
> 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, signif(mysum$coefficients[i,1],6), 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.row.start(a)
> a<-table.element(a, mywarning)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/home/pw/wessanet/rcomp/tmp/11tu551582831289.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'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,formatC(signif(mysum$coefficients[i,1],5),format='g',flag='+'))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,2],5),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,3],4),format='e',flag='+'))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,4],4),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$coefficients[i,4]/2,4),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/home/pw/wessanet/rcomp/tmp/121oa71582831289.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,formatC(signif(sqrt(mysum$r.squared),6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$r.squared,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$adj.r.squared,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a,formatC(signif(mysum$fstatistic[1],6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[2],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, signif(mysum$fstatistic[3],6))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a,formatC(signif(1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]),6),format='g',flag=' '))
> 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,formatC(signif(mysum$sigma,6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a,formatC(signif(sum(myerror*myerror),6),format='g',flag=' '))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/home/pw/wessanet/rcomp/tmp/13cf4e1582831289.tab")
> myr <- as.numeric(mysum$resid)
> myr
[1] -0.035103287 -0.051686956 0.052163892 -0.100631176 -0.191047102
[6] -0.018705550 0.013984022 -0.061270626 -0.049183723 -0.150525275
[11] -0.198022446 -0.214190190 -0.096867232 -0.154146536 -0.095724143
[16] -0.030693019 -0.111792453 0.068170360 0.174151767 0.116903181
[21] 0.116493319 0.130710780 0.145456132 0.044009083 0.235996553
[26] 0.060878120 -0.148935542 -0.107898355 -0.089389059 -0.077090756
[31] 0.026636002 0.095748372 -0.070835703 0.076859932 -0.112513149
[36] 0.087133173 -0.057773860 -0.133923011 -0.003382591 -0.085954135
[41] -0.005059632 -0.009804984 -0.015307813 -0.046637239 -0.036239907
[46] -0.074258501 0.053760093 -0.049121880 0.003636002 -0.074482025
[51] -0.082060037 -0.011296092 -0.030090756 0.077754030 -0.046389059
[56] 0.006368419 0.138231798 0.150983617 0.185332042 0.151952898
[61] -0.052407652 0.050548694 0.156157020 0.207672784 0.183548694
[66] 0.169200674
> a <-table.start()
> a <- table.row.start(a)
> a <- table.element(a,'Menu of Residual Diagnostics',2,TRUE)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <- table.element(a,'Description',1,TRUE)
> a <- table.element(a,'Link',1,TRUE)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Histogram',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_histogram.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Central Tendency',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_centraltendency.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'QQ Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_fitdistrnorm.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Kernel Density Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_density.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Skewness/Kurtosis Test',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_skewness_kurtosis.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Skewness-Kurtosis Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_skewness_kurtosis_plot.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Harrell-Davis Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_harrell_davis.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Bootstrap Plot -- Central Tendency',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_bootstrapplot1.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Blocked Bootstrap Plot -- Central Tendency',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_bootstrapplot.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'(Partial) Autocorrelation Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_autocorrelation.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Spectral Analysis',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_spectrum.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Tukey lambda PPCC Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_tukeylambda.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <-table.element(a,'Box-Cox Normality Plot',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_boxcoxnorm.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a <- table.row.start(a)
> a <- table.element(a,'Summary Statistics',1,header=TRUE)
> a <- table.element(a,hyperlink( paste('https://supernova.wessa.net/rwasp_summary1.wasp?convertgetintopost=1&data=',paste(as.character(mysum$resid),sep='',collapse=' '),sep='') ,'Compute','Click here to examine the Residuals.'),1)
> a <- table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/home/pw/wessanet/rcomp/tmp/14gj2u1582831289.tab")
> if(n < 200) {
+ 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,formatC(signif(x[i],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(x[i]-mysum$resid[i],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(mysum$resid[i],6),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/home/pw/wessanet/rcomp/tmp/15ukic1582831289.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,formatC(signif(gqarr[mypoint-kp3+1,1],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(gqarr[mypoint-kp3+1,2],6),format='g',flag=' '))
+ a<-table.element(a,formatC(signif(gqarr[mypoint-kp3+1,3],6),format='g',flag=' '))
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/home/pw/wessanet/rcomp/tmp/166rf51582831289.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,signif(numsignificant1,6))
+ a<-table.element(a,formatC(signif(numsignificant1/numgqtests,6),format='g',flag=' '))
+ 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,signif(numsignificant5,6))
+ a<-table.element(a,signif(numsignificant5/numgqtests,6))
+ 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,signif(numsignificant10,6))
+ a<-table.element(a,signif(numsignificant10/numgqtests,6))
+ 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="/home/pw/wessanet/rcomp/tmp/17y4m01582831289.tab")
+ }
+ }
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Ramsey RESET F-Test for powers (2 and 3) of fitted values',1,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> reset_test_fitted <- resettest(mylm,power=2:3,type='fitted')
> a<-table.element(a,paste('
',RC.texteval('reset_test_fitted'),'',sep='')) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Ramsey RESET F-Test for powers (2 and 3) of regressors',1,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > reset_test_regressors <- resettest(mylm,power=2:3,type='regressor') > a<-table.element(a,paste('
',RC.texteval('reset_test_regressors'),'',sep='')) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'Ramsey RESET F-Test for powers (2 and 3) of principal components',1,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > reset_test_principal_components <- resettest(mylm,power=2:3,type='princomp') > a<-table.element(a,paste('
',RC.texteval('reset_test_principal_components'),'',sep='')) > a<-table.row.end(a) > a<-table.end(a) > table.save(a,file="/home/pw/wessanet/rcomp/tmp/18bb4t1582831289.tab") > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Variance Inflation Factors (Multicollinearity)',1,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > vif <- vif(mylm) Error in vif.default(mylm) : model contains fewer than 2 terms Calls: vif -> vif.default Execution halted