R Under development (unstable) (2024-11-24 r87369 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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.

> 
> #""""""""""""""""""""""""Detecting periodicity
> #install.packages("expm", repos = "https://cran.r-project.org")
> #install.packages("sn", repos = "https://cran.r-project.org")
> #install.packages("readxl", repos = "https://cran.r-project.org")
> 
> #install.packages("expm")
> #install.packages("sn")
> #install.packages("readxl")
> 
> 
> #library(PerRegMod)
> requireNamespace("PerRegMod")
Loading required namespace: PerRegMod
> requireNamespace("expm")
> requireNamespace("readxl")
> requireNamespace("sn")
> 
> library(PerRegMod)
> library(expm)
Loading required package: Matrix

Attaching package: 'expm'

The following object is masked from 'package:Matrix':

    expm

> library(readxl)
> library(sn)# for skew normal density
Loading required package: stats4

Attaching package: 'sn'

The following object is masked from 'package:stats':

    sd

> 
> 
> 
> n=200
> s=2
> x1=rnorm(n,0,1)
> x2=rnorm(n,0,2)
> x3=rnorm(n,0,3)
> x4=rnorm(n,0,2.7)
> y=rnorm(n,0,2.5)
> x=list(x1,x2,x3,x4)
> model=lm(y~x1+x2+x3+x4)
> z=model$residuals
> check_periodicity(x,y,s)
Chi-squared test for detecting periodicity of coefficients
Chi-squared statistic: 27.3695  on df= 6 , p-value:  0.0001235> 
> ####""""""""""""""""""""""""" LSE method
> set.seed(4)
> s=2
> n=200
> m=n/s
> p=3
> mu=c(2,6)
> beta1=c(3,5)
> beta2=c(1,2.5)
> beta3=c(-2,1)
> x1=runif(n,0,5)
> x2=runif(n,0,10)
> x3=runif(n,0,15)
> y=rep(0,n)
> for (i in 1:s) {
+   q=seq(i,n,s)
+   y[q]=mu[i] + beta1[i] * x1[q] + beta2[i] * x2[q] + beta3[i] * x3[q] + rsn(m,alpha=10)
+ }
> x=list(x1,x2,x3)
> lm_per(x,y,s)
Residuals:
      0%      25%      50%      75%     100% 
-0.93521 -0.45666 -0.13735  0.33775  2.19037 

Coefficients:
     intercept x...1... X...2... X...3...      sd
s= 1    2.8421   2.9729  0.99535  -1.9957 0.61037
s= 2    6.7879   4.9639  2.51007   1.0028 0.57218

R-squared:  0.997132 Adjusted R-squared:  0.997042
Root mean square error:  0.5796334> 
> 
> ##""""""""""""""""""""""" AE
> set.seed(2)
> s=2
> n=200
> m=n/s
> p=3
> mu=c(2,6)
> beta1=c(3,5)
> beta2=c(1,2.5)
> beta3=c(-2,1)
> x1=runif(n,0,5)
> x2=runif(n,0,10)
> x3=runif(n,0,15)
> y=rep(0,n)
> for (i in 1:s) {
+   q=seq(i,n,s)
+   y[q]=mu[i] + beta1[i] * x1[q] + beta2[i] * x2[q] + beta3[i] * x3[q] + rnorm(m,0,1)
+ }
> x=list(x1,x2,x3)
> lm_per_AE(x,y,s)
Residuals:
         0%         25%         50%         75%        100% 
-2.58606414 -0.63802606  0.00068852  0.65515144  2.11237931 

Coefficients:
     intercept x...1... X...2... X...3...      SD
s= 1    2.0033   3.0217   1.0090  -1.9751 0.80586
s= 2    6.0168   5.0314   2.4428   1.0017 0.98376

Root mean square error:  0.9355969> 
> ######"""""""""""""""""""""""""" real data examples
> #####"
> #Data <- system.file("data", "weather_data_casablanca-settat.xlsx", package = "PerRegMod")
> #load("~/PerRegMod_4.4.1/PerRegMod/data/weather_data_casablanca_settat.RData")
> #Data=weather_data_casablanca_settat
> 
> #file_path <- system.file("data", "weather_data_casablanca_settat.xlsx", package = "PerRegMod")
> #Data <- readxl::read_excel(file_path)
> #file_path <- system.file("inst/extdata", "weather_data_casablanca_settat.xlsx", package = "PerRegMod")
> #Data <- readxl::read_excel(file_path)
> #Data=read_xlsx("inst/extdata/weather_data_casablanca_settat.xlsx")
> file_path <- system.file("extdata", "weather_data_casablanca_settat.xlsx", package = "PerRegMod")
> 
> # Check if the file path is found, then read the file
> if (file_path != "") {
+   Data <- readxl::read_excel(file_path)
+ } else {
+   stop("File not found. Please check that the file is correctly placed in 'inst/extdata' within the package directory.")
+ }
> months <- c('January','February','March','April','May','June',
+             'July','August','September','October','November','December')
> years <- c('2010','2011','2012','2013','2014','2015','2016',
+            '2017','2018','2019','2020')
> 
> #########""""""""""""for temperature
> y=rep(0,11 * 12)
> ly=length(years)
> lm=length(months)
> for (i in 1:ly ) {
+   for (j in 1:lm) {
+     y[(i - 1) * lm + j] = mean(Data$Temperature..C[Data$Year == years[i] &
+                                                      Data$Month == months[j] ]  )
+   }
+ 
+ }
> 
> #########""""""""""""for Humidity
> x1=rep(0,11 * 12)
> ly=length(years)
> lm=length(months)
> for (i in 1:ly ) {
+   for (j in 1:lm) {
+     x1[(i - 1) * lm + j]=mean(Data$Humidity..[Data$Year==years[i] &
+                                                 Data$Month==months[j] ]  )
+   }
+ 
+ }
> 
> #########""""""""""""for dew point
> x2=rep(0,11 * 12)
> ly=length(years)
> lm=length(months)
> for (i in 1:ly ) {
+   for (j in 1:lm) {
+     x2[(i - 1) * lm + j]=mean(Data$Dew.Point..C[Data$Year==years[i] &
+                                                   Data$Month==months[j] ]  )
+   }
+ 
+ }
> 
> #########""""""""""""for wind
> x3=rep(0,11 * 12)
> ly=length(years)
> lm=length(months)
> for (i in 1:ly ) {
+   for (j in 1:lm) {
+     x3[(i - 1) * lm + j]=mean(Data$Wind.Speed.km.h[Data$Year==years[i] &
+                                                      Data$Month==months[j] ]  )
+   }
+ 
+ }
> 
> ###""""""""""""""""""""""""example 1
> 
> #""""""""""check periodicity
> x=list(x2,x3)# dew point and wind
> n=length(y)
> s=12
> m=n/s
> check_periodicity(x,y,s)
Chi-squared test for detecting periodicity of coefficients
Chi-squared statistic: 139.35296  on df= 44 , p-value:  0> 
> #""""""""estimating parameters using LSE
> lm_per(x,y,s)
Residuals:
       0%       25%       50%       75%      100% 
-3.797334 -0.725823 -0.023063  0.731557  3.387731 

Coefficients:
      intercept  x...1...   X...2...      sd
s= 1    14.7686  0.365816 -0.3447470 0.77248
s= 2    14.4031  0.783058 -0.6049130 1.31027
s= 3    20.2313 -0.095313 -0.2104771 1.01714
s= 4    16.6927  0.893365 -0.6399310 1.44247
s= 5    13.5495  0.767417  0.0056104 2.07729
s= 6     5.1898  1.600426 -0.3040881 1.06302
s= 7    42.1296 -0.431045 -0.6816115 1.68553
s= 8    37.8579 -0.247980 -0.4538573 1.31940
s= 9    17.0024  0.437180  0.2094528 0.77024
s= 10   32.7012 -0.306806 -0.5037246 1.72124
s= 11   26.9001 -0.489266 -0.2760296 1.26232
s= 12   19.4080 -0.071502 -0.2781806 1.15929

R-squared:  0.341975 Adjusted R-squared:  0.177469
Root mean square error:  1.153366> #""""""""estimating parameters using AE
> lm_per_AE(x,y,s)
Residuals:
      0%      25%      50%      75%     100% 
-3.91598 -0.71694 -0.02005  0.69440  3.40537 

Coefficients:
      intercept  x...1...  X...2...      SD
s= 1    14.7685  0.365761 -0.344803 0.59670
s= 2    14.4032  0.783071 -0.604987 1.72700
s= 3    20.2313 -0.095259 -0.210404 1.03221
s= 4    16.6924  0.893116 -0.640378 2.08507
s= 5    13.5540  0.771908  0.010586 4.38878
s= 6     5.1896  1.600153 -0.304393 1.13234
s= 7    42.1290 -0.431800 -0.682552 2.83709
s= 8    37.8577 -0.248023 -0.453798 1.71883
s= 9    17.0026  0.437445  0.209642 0.59308
s= 10   32.7003 -0.307424 -0.504580 2.92165
s= 11   26.8998 -0.489503 -0.276246 1.59654
s= 12   19.4079 -0.071575 -0.278251 1.34447

Root mean square error:  1.153864> 
> ####"""""""""""""""""""""Example 2
> ##""""""check periodicity
> x=list(x1,x3)# humidity and wind
> n=length(y)
> s=12
> m=n/s
> check_periodicity(x,y,s)
Chi-squared test for detecting periodicity of coefficients
Chi-squared statistic: 1131.78449  on df= 44 , p-value:  0> 
> #""""""""estimating parameters using LSE
> lm_per(x,y,s)
Residuals:
       0%       25%       50%       75%      100% 
-2.368374 -0.479260 -0.083922  0.544323  3.561058 

Coefficients:
      intercept  x...1...  X...2...      sd
s= 1     14.980  0.033304 -0.313531 0.93041
s= 2     20.710 -0.089701  0.122883 1.90827
s= 3     29.059 -0.156965 -0.061194 0.64481
s= 4     46.467 -0.354884 -0.268507 0.82568
s= 5     43.157 -0.306016 -0.115093 1.51031
s= 6     50.872 -0.378169 -0.225892 0.99913
s= 7     48.472 -0.322814 -0.140465 1.00970
s= 8     45.152 -0.250712 -0.179926 0.68560
s= 9     44.828 -0.262108 -0.236498 0.57371
s= 10    42.237 -0.232873 -0.458102 0.63577
s= 11    32.466 -0.163448 -0.290638 0.75271
s= 12    24.966 -0.106498 -0.155004 0.96201

R-squared:  0.622292 Adjusted R-squared:  0.527865
Root mean square error:  0.8738242> #""""""""estimating parameters using AE
> lm_per_AE(x,y,s)
Residuals:
      0%      25%      50%      75%     100% 
-2.50805 -0.48476  0.01095  0.59413  4.38712 

Coefficients:
      intercept  x...1...  X...2...       SD
s= 1     14.946  0.032701 -0.314299  0.76555
s= 2     20.067 -0.092184  0.121675 -4.63692
s= 3     29.084 -0.156333 -0.060532  0.45120
s= 4     46.464 -0.354666 -0.268246  0.67259
s= 5     43.344 -0.307257 -0.112224  3.79255
s= 6     51.002 -0.377706 -0.224841  1.46153
s= 7     48.405 -0.324608 -0.141179  0.78229
s= 8     45.153 -0.251054 -0.180911  0.47392
s= 9     44.834 -0.260922 -0.235151  0.33383
s= 10    42.198 -0.234702 -0.459730  0.35070
s= 11    32.460 -0.161397 -0.288685  0.54968
s= 12    24.990 -0.104948 -0.153958  0.99888

Root mean square error:  0.9143536> 
> 
> proc.time()
   user  system elapsed 
   4.40    0.21    4.62