unit root test results output

[kdom20]

Hi all,

1st time using R so please bear with me.
nonlinearTSA is a new package to run some nonlinear unit root tests, i.e., Sollis 2004.
the function runs fine and is:
Sollis_2004_unit_root(x, model = 1, max_lags = 6)
Since I have many time series (169 in total) and to save time, I utilized the apply function as below and again it runs fine.

apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)

the two questions I have are:

1. the output is not printing the series' names!
2. how can I export the output/results into excel.

many thanks.
Karla

[HanOostdijk]

Below I give an example of what you could do.
Not having your data set I create my own consisting of num_series (2) time series.
The time series are named c1 and c2.
Then I :

• use your code and assign the result to variable y . y is a named list with the model results.
• would like to apply on each model result the glance and tidy functions of the broom package but this model type is apparently not supported (yet).
• print the model result with the name of model (column name). This answers your first question.
• show how to extract and print the coefficients (**) of one model. I have to code this function myself because broom does not support these models.
• expand this function to extract the coefficients (**) of all models into a data.frame .
• write this data.frame to a csv file so that the data can be included in an Excel file

(**) Maybe you want also to copy other things (e.g. the residuals). In that case you will have to analyze the contents of an y element and when you have found the item copy it to a data.frame like done here for the coefficients. As a start for the analysis :

names(y[[1]]$Model)
[1] "object" "mse" "AIC" "MAPE" "df" "residuals"
[7] "lowCoef" "highCoef" "thCoef" "fixedTh" "externThVar" "coef"
[13] "mTh"

Code to extract and save coefficients in a csv file:

library(NonlinearTSA)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)
set.seed(123)

num_series <- 2
x <-  rnorm(1000*num_series)
da <-  matrix(x,nrow=1000,ncol=num_series,byrow = F,
            dimnames=list(NULL,paste('c',1:num_series,sep = '')))
str(da)
#>  num [1:1000, 1:2] -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : chr [1:2] "c1" "c2"

y  <-  apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
names(y) 
#> [1] "c1" "c2"

# broom::glance(y[[1]]) # no broom routine available (?)

# print model results with the name of model
purrr::iwalk(y, function(yd,yn) {
      cat(paste('\n Result model',yn,":\n\n" ))
      print(yd)
    }
  )
#> 
#>  Result model c1 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>      phiL.1      phiH.1     Dphi. 1 
#> -1.05658408 -1.06565697  0.03141912 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 50.2%     High regime: 49.8% 
#> 
#> Residuals:
#>        Min         1Q     Median         3Q        Max 
#> -2.7813568 -0.6500425 -0.0017678  0.6508761  3.1574841 
#> 
#> Fit:
#> residuals variance = 0.979,  AIC = -15, MAPE = 733.2%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value Pr(>|t|)    
#> phiL.1  -1.056584    0.055395 -19.0737   <2e-16 ***
#> phiH.1  -1.065657    0.055306 -19.2683   <2e-16 ***
#> Dphi. 1  0.031419    0.031658   0.9925   0.3212    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> $`Selected lag`
#> [1] 1
#> 
#> $`p1=p2=0 Statistic`
#> [1] 273.0318
#> 
#> 
#>  Result model c2 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>      phiL.1      phiH.1     Dphi. 1 
#> -1.07626502 -1.08815727  0.03779523 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 49.7%     High regime: 50.3% 
#> 
#> Residuals:
#>        Min         1Q     Median         3Q        Max 
#> -3.1398771 -0.6964508  0.0043722  0.7196956  3.3904897 
#> 
#> Fit:
#> residuals variance = 1.012,  AIC = 18, MAPE = 254.7%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value Pr(>|t|)    
#> phiL.1  -1.076265    0.055447 -19.4108   <2e-16 ***
#> phiH.1  -1.088157    0.055585 -19.5763   <2e-16 ***
#> Dphi. 1  0.037795    0.031618   1.1954   0.2322    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> $`Selected lag`
#> [1] 1
#> 
#> $`p1=p2=0 Statistic`
#> [1] 281.1227

# example for extracting (only) the coefficients for one (first) model 
df1  <-  as.data.frame(y[[1]]$Model$coef) %>%
          mutate(var=rownames(.))
df1  <- cbind(data.frame(modelname=names(y)[1]),df1 ) 
print(df1)
#>   modelname    Estimate  Std. Error     t value     Pr(>|t|)     var
#> 1        c1 -1.05658408  0.05539470 -19.0737393 2.185959e-69  phiL.1
#> 2        c1 -1.06565697  0.05530630 -19.2682733 1.422443e-70  phiH.1
#> 3        c1  0.03141912  0.03165804   0.9924532 3.212173e-01 Dphi. 1

# example for extracting (only) the coefficients for all models
df_all <- purrr::imap_dfr(y, function(yd,yn) {
      df  <-  as.data.frame(yd$Model$coef) %>%
          mutate(var=rownames(.))
      cbind(data.frame(modelname=yn),df )
    }
  )

print(df_all)
#>   modelname    Estimate  Std. Error     t value     Pr(>|t|)     var
#> 1        c1 -1.05658408  0.05539470 -19.0737393 2.185959e-69  phiL.1
#> 2        c1 -1.06565697  0.05530630 -19.2682733 1.422443e-70  phiH.1
#> 3        c1  0.03141912  0.03165804   0.9924532 3.212173e-01 Dphi. 1
#> 4        c2 -1.07626502  0.05544685 -19.4107528 1.907635e-71  phiL.1
#> 5        c2 -1.08815727  0.05558530 -19.5763485 1.831471e-72  phiH.1
#> 6        c2  0.03779523  0.03161848   1.1953527 2.322335e-01 Dphi. 1
# write coefficients for all models to csv file (to be read in xls file)
write.csv(df_all,'output.csv') # see ?write.csv for more arguments if necessary

^{Created on 2020-07-18 by the reprex package (v0.3.0)}

[kdom20]

@HanOostdijk

library(NonlinearTSA)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)
set.seed(123)

num_series <- 2
x <-  rnorm(1000*num_series)
da <-  matrix(x,nrow=1000,ncol=num_series,byrow = F,
            dimnames=list(NULL,paste('c',1:num_series,sep = '')))
str(da)
#>  num [1:1000, 1:2] -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : chr [1:2] "c1" "c2"

y  <-  apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
names(y) 
#> [1] "c1" "c2"

# broom::glance(y[[1]]) # no broom routine available (?)

# print model results with the name of model
purrr::iwalk(y, function(yd,yn) {
      cat(paste('\n Result model',yn,":\n\n" ))
      print(yd)
    }
  )
#> 
#>  Result model c1 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>      phiL.1      phiH.1     Dphi. 1 
#> -1.05658408 -1.06565697  0.03141912 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 50.2%     High regime: 49.8% 
#> 
#> Residuals:
#>        Min         1Q     Median         3Q        Max 
#> -2.7813568 -0.6500425 -0.0017678  0.6508761  3.1574841 
#> 
#> Fit:
#> residuals variance = 0.979,  AIC = -15, MAPE = 733.2%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value Pr(>|t|)    
#> phiL.1  -1.056584    0.055395 -19.0737   <2e-16 ***
#> phiH.1  -1.065657    0.055306 -19.2683   <2e-16 ***
#> Dphi. 1  0.031419    0.031658   0.9925   0.3212    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> ```
Selected lag`
#> [1] 1
#> 
#> ```
p1=p2=0 Statistic`
#> [1] 273.0318
#> 
#> 
#>  Result model c2 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>      phiL.1      phiH.1     Dphi. 1 
#> -1.07626502 -1.08815727  0.03779523 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 49.7%     High regime: 50.3% 
#> 
#> Residuals:
#>        Min         1Q     Median         3Q        Max 
#> -3.1398771 -0.6964508  0.0043722  0.7196956  3.3904897 
#> 
#> Fit:
#> residuals variance = 1.012,  AIC = 18, MAPE = 254.7%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value Pr(>|t|)    
#> phiL.1  -1.076265    0.055447 -19.4108   <2e-16 ***
#> phiH.1  -1.088157    0.055585 -19.5763   <2e-16 ***
#> Dphi. 1  0.037795    0.031618   1.1954   0.2322    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> ```
Selected lag`
#> [1] 1
#> 
#> ```
p1=p2=0 Statistic`
#> [1] 281.1227

# example for extracting (only) the coefficients for one (first) model 
df1  <-  as.data.frame(y[[1]]$Model$coef) %>%
          mutate(var=rownames(.))
df1  <- cbind(data.frame(modelname=names(y)[1]),df1 ) 
print(df1)
#>   modelname    Estimate  Std. Error     t value     Pr(>|t|)     var
#> 1        c1 -1.05658408  0.05539470 -19.0737393 2.185959e-69  phiL.1
#> 2        c1 -1.06565697  0.05530630 -19.2682733 1.422443e-70  phiH.1
#> 3        c1  0.03141912  0.03165804   0.9924532 3.212173e-01 Dphi. 1

# example for extracting (only) the coefficients for all models
df_all <- purrr::imap_dfr(y, function(yd,yn) {
      df  <-  as.data.frame(yd$Model$coef) %>%
          mutate(var=rownames(.))
      cbind(data.frame(modelname=yn),df )
    }
  )

print(df_all)
#>   modelname    Estimate  Std. Error     t value     Pr(>|t|)     var
#> 1        c1 -1.05658408  0.05539470 -19.0737393 2.185959e-69  phiL.1
#> 2        c1 -1.06565697  0.05530630 -19.2682733 1.422443e-70  phiH.1
#> 3        c1  0.03141912  0.03165804   0.9924532 3.212173e-01 Dphi. 1
#> 4        c2 -1.07626502  0.05544685 -19.4107528 1.907635e-71  phiL.1
#> 5        c2 -1.08815727  0.05558530 -19.5763485 1.831471e-72  phiH.1
#> 6        c2  0.03779523  0.03161848   1.1953527 2.322335e-01 Dphi. 1
# write coefficients for all models to csv file (to be read in xls file)
write.csv(df_all,'output.csv') # see ?write.csv for more arguments if necessary

HanOostddijk,
This is incredible and I can't thank you enough for taking the time. I will run your code on my data asap and if I have any more questions, I hope you do not mind me following up.
Again, many thanks.
K

[kdom20]

library(NonlinearTSA)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)

num_series <- 2
da <- GDP1_2
#> Error in eval(expr, envir, enclos): object 'GDP1_2' not found
da <-  matrix(x,nrow=317,ncol=num_series,byrow = F,
              dimnames=list(NULL,paste('c',1:num_series,sep = '')))
#> Error in matrix(x, nrow = 317, ncol = num_series, byrow = F, dimnames = list(NULL, : object 'x' not found

y  <-  apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
#> Error in apply(da, 2, Sollis_2004_unit_root, model = 3, max_lags = 6): object 'da' not found
names(y) 
#> Error in eval(expr, envir, enclos): object 'y' not found
str(da)
#> Error in str(da): object 'da' not found
purrr::iwalk(y, function(yd,yn) {
    cat(paste('\n Result model',yn,":\n\n" ))
    print(yd)
}
)
#> Error in map2(.x, .y, .f, ...): object 'y' not found

df1  <-  as.data.frame(y[[1]]$Model$coef) %>%
    mutate(var=rownames(.))
#> Error in as.data.frame(y[[1]]$Model$coef): object 'y' not found
df1  <- cbind(data.frame(modelname=names(y)[1]),df1 ) 
#> Error in data.frame(modelname = names(y)[1]): object 'y' not found
print(df1)
#> Error in print(df1): object 'df1' not found


df2  <-  as.data.frame(y[[1]]$Model$coef) %>%
    mutate(var=rownames(.))
#> Error in as.data.frame(y[[1]]$Model$coef): object 'y' not found
df1  <- cbind(data.frame(modelname=names(y)[1]),df1 ) 
#> Error in data.frame(modelname = names(y)[1]): object 'y' not found
print(df1)
#> Error in print(df1): object 'df1' not found



df_all <- purrr::imap_dfr(y, function(yd,yn) {
    df  <-  as.data.frame(yd$Model$coef) %>%
        mutate(var=rownames(.))
    cbind(data.frame(modelname=yn),df )
}
)
#> Error in map2(.x, .y, .f, ...): object 'y' not found

print(df_all)
#> Error in print(df_all): object 'df_all' not found

write.csv(df_all,"C:\\ENERGY_169\\Results\\output1.csv", row.names = T)
#> Error in is.data.frame(x): object 'df_all' not found

^{Created on 2020-07-18 by the reprex package (v0.3.0)}

[kdom20]

Trying to use reprex but obviously it did not work as intended. My apology.
The code seems to be running fine but with the error! (especially the two persistent ones, i.e., object 'x' not found and object 'res' not found.

I tried to upload couple of series of the data but it seems all types of data files are not allowed to be uploaded!

[kdom20]

not sure if this will work at your end, but this is a paste option for couple of series with just T = 100. There are other pate options but I have no clue which works and which does not. I also tried to create an additional DF to extract just the selected lag and the p1=p2=0 test Statistic and that failed miserably!
If this will require more of your time, I totally understand a no response. Thanks again.

data.frame(
stringsAsFactors = FALSE,
series1 = c("series2","1","10364.956",
"9036.963","2","11294.906","9847.764","3","11003.601",
"9593.782","4","10298.415","8978.948","5",
"12260.282","10689.453","6","11636.513","10145.603","7",
"9854.315","8591.746","8","11071.253","9652.766","9",
"11784.406","10274.548","10","10209.250","8901.207",
"11","9281.871","8092.646","12","10120.241","8823.601",
"13","10110.924","8815.478","14","9676.897",
"8437.060","15","10644.657","9280.827","16","10345.933",
"9020.377","17","10737.304","9361.604","18",
"11342.680","9889.417","19","11983.218","10447.888","20",
"11772.107","10263.825","21","12522.962","10918.478",
"22","12034.415","10492.525","23","12162.756",
"10604.422","24","11870.582","10349.683","25","11752.067",
"10246.352","26","12058.674","10513.675","27",
"11686.437","10189.131","28","11615.530","10127.309","29",
"11969.968","10436.336","30","11304.754","9856.350",
"31","11563.501","10081.946","32","11379.801",
"9921.782","33","12135.886","10580.995","34","12404.936",
"10815.573","35","12360.136","10776.514","36",
"12066.902","10520.850","37","12963.839","11302.868","38",
"12269.519","10697.506","39","12498.734","10897.354",
"40","12415.221","10824.541","41","12206.651",
"10642.693","42","12563.076","10953.452","43",
"12852.312","11205.631","44","12550.040","10942.086","45",
"12374.617","10789.139","46","12556.925","10948.089",
"47","12951.392","11292.016","48","13541.552",
"11806.563","49","13396.669","11680.243","50","13317.196",
"11610.952","51","13510.226","11779.250","52",
"13217.920","11524.395","53","13868.934","12091.999","54",
"13992.017","12199.312","55","13599.989","11857.512",
"56","14012.198","12216.907","57","13774.808",
"12009.933","58","14568.218","12701.689","59","14770.864",
"12878.371","60","15004.525","13082.094","61",
"15440.562","13462.264","62","16055.198","13998.151",
"63","15608.598","13608.771","64","15557.349",
"13564.089","65","15550.734","13558.321","66","15424.645",
"13448.387","67","15630.701","13628.043","68",
"15738.359","13721.907","69","15883.316","13848.292","70",
"16780.869","14630.847","71","16322.372","14231.095",
"72","16627.296","14496.951","73","16202.553",
"14126.627","74","16491.447","14378.507","75","16124.602",
"14058.664","76","16639.779","14507.834","77",
"17178.750","14977.750","78","17675.692","15411.022","79",
"17552.409","15303.535","80","17504.365",
"15261.647","81","18304.071","15958.892","82","18825.509",
"16413.521","83","18637.946","16249.989","84",
"18542.012","16166.346","85","18615.162","16230.124","86",
"18485.428","16117.012","87","18808.771","16398.927",
"88","18495.652","16125.926","89","18453.289",
"16088.991","90","18866.279","16449.067","91","20298.118",
"17697.454","92","20385.070","17773.265","93",
"21204.286","18487.521","94","20860.450","18187.738","95",
"20185.320","17599.108","96","22377.818","19510.696",
"97","22779.405","19860.831","98","22614.071",
"19716.679","99","22831.421","19906.182","100",
"23431.460","20429.341")
)

[HanOostdijk]

In a previous entry the error was that the (data.frame ?) GDP1_2 could not be found and a line later you try to use x that was not defined.
In the last entry the format of your da data.frame was not good: I corrected it in the code.

library(NonlinearTSA)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)
set.seed(123)

da = data.frame(
stringsAsFactors = FALSE,
series1 = c("series2","1","10364.956",
"9036.963","2","11294.906","9847.764","3","11003.601",
"9593.782","4","10298.415","8978.948","5",
"12260.282","10689.453","6","11636.513","10145.603","7",
"9854.315","8591.746","8","11071.253","9652.766","9",
"11784.406","10274.548","10","10209.250","8901.207",
"11","9281.871","8092.646","12","10120.241","8823.601",
"13","10110.924","8815.478","14","9676.897",
"8437.060","15","10644.657","9280.827","16","10345.933",
"9020.377","17","10737.304","9361.604","18",
"11342.680","9889.417","19","11983.218","10447.888","20",
"11772.107","10263.825","21","12522.962","10918.478",
"22","12034.415","10492.525","23","12162.756",
"10604.422","24","11870.582","10349.683","25","11752.067",
"10246.352","26","12058.674","10513.675","27",
"11686.437","10189.131","28","11615.530","10127.309","29",
"11969.968","10436.336","30","11304.754","9856.350",
"31","11563.501","10081.946","32","11379.801",
"9921.782","33","12135.886","10580.995","34","12404.936",
"10815.573","35","12360.136","10776.514","36",
"12066.902","10520.850","37","12963.839","11302.868","38",
"12269.519","10697.506","39","12498.734","10897.354",
"40","12415.221","10824.541","41","12206.651",
"10642.693","42","12563.076","10953.452","43",
"12852.312","11205.631","44","12550.040","10942.086","45",
"12374.617","10789.139","46","12556.925","10948.089",
"47","12951.392","11292.016","48","13541.552",
"11806.563","49","13396.669","11680.243","50","13317.196",
"11610.952","51","13510.226","11779.250","52",
"13217.920","11524.395","53","13868.934","12091.999","54",
"13992.017","12199.312","55","13599.989","11857.512",
"56","14012.198","12216.907","57","13774.808",
"12009.933","58","14568.218","12701.689","59","14770.864",
"12878.371","60","15004.525","13082.094","61",
"15440.562","13462.264","62","16055.198","13998.151",
"63","15608.598","13608.771","64","15557.349",
"13564.089","65","15550.734","13558.321","66","15424.645",
"13448.387","67","15630.701","13628.043","68",
"15738.359","13721.907","69","15883.316","13848.292","70",
"16780.869","14630.847","71","16322.372","14231.095",
"72","16627.296","14496.951","73","16202.553",
"14126.627","74","16491.447","14378.507","75","16124.602",
"14058.664","76","16639.779","14507.834","77",
"17178.750","14977.750","78","17675.692","15411.022","79",
"17552.409","15303.535","80","17504.365",
"15261.647","81","18304.071","15958.892","82","18825.509",
"16413.521","83","18637.946","16249.989","84",
"18542.012","16166.346","85","18615.162","16230.124","86",
"18485.428","16117.012","87","18808.771","16398.927",
"88","18495.652","16125.926","89","18453.289",
"16088.991","90","18866.279","16449.067","91","20298.118",
"17697.454","92","20385.070","17773.265","93",
"21204.286","18487.521","94","20860.450","18187.738","95",
"20185.320","17599.108","96","22377.818","19510.696",
"97","22779.405","19860.831","98","22614.071",
"19716.679","99","22831.421","19906.182","100",
"23431.460","20429.341")
)


# da has not the correct format (look at it) therefore correct
series1 = as.numeric(da$series1[seq(3,length(da$series1),3)])
series2 = as.numeric(da$series1[seq(4,length(da$series1),3)])
da_corrected = data.frame(series1=series1,series2=series2)
# this is the correct format (look at difference with da)

y  <-  apply(da_corrected,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
names(y) 
#> [1] "series1" "series2"

# broom::glance(y[[1]]) # no broom routine available (?)

# print model results with the name of model
purrr::iwalk(y, function(yd,yn) {
      cat(paste('\n Result model',yn,":\n\n" ))
      print(yd)
    }
  )
#> 
#>  Result model series1 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>     phiL.1     phiH.1    Dphi. 1 
#> -0.5009353 -0.5353598  0.0110481 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 46.94%    High regime: 53.06% 
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -1307.019  -275.034    19.989   277.564  1775.622 
#> 
#> Fit:
#> residuals variance = 248558,  AIC = 1248, MAPE = 114.6%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value  Pr(>|t|)    
#> phiL.1  -0.500935    0.130864  -3.8279 0.0002290 ***
#> phiH.1  -0.535360    0.139976  -3.8246 0.0002316 ***
#> Dphi. 1  0.011048    0.100733   0.1097 0.9128924    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> $`Selected lag`
#> [1] 1
#> 
#> $`p1=p2=0 Statistic`
#> [1] 12.73111
#> 
#> 
#>  Result model series2 :
#> 
#> $Model
#> 
#> Non linear autoregressive model
#> 
#> SETAR model ( 2 regimes)
#> Coefficients:
#>      phiL.1      phiH.1     Dphi. 1 
#> -0.50116068 -0.53675454  0.01133583 
#> 
#> Threshold:
#> -Variable: Z(t) = + (1) X(t)
#> -Value: 0 (fixed)
#> Proportion of points in low regime: 46.94%    High regime: 53.06% 
#> 
#> Residuals:
#>       Min        1Q    Median        3Q       Max 
#> -1137.958  -239.844    16.838   242.025  1548.452 
#> 
#> Fit:
#> residuals variance = 188820,  AIC = 1221, MAPE = 114.8%
#> 
#> Coefficient(s):
#> 
#>          Estimate  Std. Error  t value  Pr(>|t|)    
#> phiL.1  -0.501161    0.130811  -3.8312 0.0002263 ***
#> phiH.1  -0.536755    0.140081  -3.8317 0.0002259 ***
#> Dphi. 1  0.011336    0.100710   0.1126 0.9106127    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Threshold
#> Variable: Z(t) = + (1) X(t) 
#> 
#> Value: 0 (fixed)
#> 
#> $`Selected lag`
#> [1] 1
#> 
#> $`p1=p2=0 Statistic`
#> [1] 12.76459

# example for extracting (only) the coefficients for one (first) model 
df1  <-  as.data.frame(y[[1]]$Model$coef) %>%
          mutate(var=rownames(.))
df1  <- cbind(data.frame(modelname=names(y)[1]),df1 ) 
print(df1)
#>   modelname   Estimate  Std. Error    t value     Pr(>|t|)     var
#> 1   series1 -0.5009353   0.1308643 -3.8278998 0.0002289815  phiL.1
#> 2   series1 -0.5353598   0.1399763 -3.8246472 0.0002316216  phiH.1
#> 3   series1  0.0110481   0.1007333  0.1096767 0.9128923554 Dphi. 1

# example for extracting (only) the coefficients for all models
df_all <- purrr::imap_dfr(y, function(yd,yn) {
      df  <-  as.data.frame(yd$Model$coef) %>%
          mutate(var=rownames(.))
      cbind(data.frame(modelname=yn),df )
    }
  )

print(df_all)
#>   modelname    Estimate  Std. Error    t value     Pr(>|t|)     var
#> 1   series1 -0.50093527   0.1308643 -3.8278998 0.0002289815  phiL.1
#> 2   series1 -0.53535980   0.1399763 -3.8246472 0.0002316216  phiH.1
#> 3   series1  0.01104810   0.1007333  0.1096767 0.9128923554 Dphi. 1
#> 4   series2 -0.50116068   0.1308108 -3.8311882 0.0002263414  phiL.1
#> 5   series2 -0.53675454   0.1400813 -3.8317367 0.0002259039  phiH.1
#> 6   series2  0.01133583   0.1007100  0.1125591 0.9106127080 Dphi. 1
# write coefficients for all models to csv file (to be read in xls file)
write.csv(df_all,'output.csv') # see ?write.csv for more arguments if necessary

^{Created on 2020-07-19 by the reprex package (v0.3.0)}

[kdom20]

Great. Many thanks. this works as well.
I think this happened because of how I pasted the data (da).
However, instead of using the data.frame for da as in your last code, will 'attach" work? and how to define x in this instance? My main data (GDP1_2) that I am reading in the program has many series (series1, series2, series3, ... etc). Something like the below just ran on just two series and it ran with no errors and that, to be honest, concerned me! as I was not sure how to define x!. Of course, I also compared the results with just running the original function on individual series and the results match!

library(NonlinearTSA)
library(dplyr)
library(purrr)

attach(GDP1_2)
num_series <- 2
da <- GDP1_2
da <- matrix(x,nrow=317,ncol=num_series,byrow = F,
dimnames=list(NULL,paste('c',1:num_series,sep = '')))

y <- apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
names(y)

I promise the above is the the last question.
K

[HanOostdijk]

Sorry for the late response and probably you have found the answer yourself.
If that is not the case:

There is no need to anything with x and also no need to 'attach' GDP1_2 .

When you have GDP1_2 in your session available as a data.frame you can run the code

da <- GDP1_2
y <- apply(da,2,Sollis_2004_unit_root,model = 3, max_lags = 6)
and/with the rest of the code

Of course that supposes that GDP1_2 only contains the time-series in numerical format as columns/fields. Otherwise the data.frame has to converted to the right format.

If GDP1_2 is not yet available as a data.frame, you have to use some code to get it in that format.
E.g. when the data is available in a xls-files, you can save it as an csv-file (e.g. GDP1_2.csv') and read the latter file with code like e.g.
GDP1_2 = read.csv('GDP1_2.csv')
If you need to read the data from a file and have no code yet to read it use the previous line.
If you have problems with that, show us the first two lines of the csv-file.
If the input is not in an excel file, let us know.

[kdom20]

Thank you so much. , I really appreciate you taking the time. As a WinRATS user, you definitely made made my first R session less stressful.

Thanks.

K