source("gh-final.R")
xi <- 0
omega <- 2
g <- 0.5
h <- 0.3

n <- 400

z <- rnorm(n)
#generate g-h data
y <- gh_trans(z,g=g,h=h)*omega+xi

nbreaks<-max(n,1000) ##number of knots


#fit the g-and-h model                    
obj <- gh_likest_linear(y,par0=NULL,nbreaks=nbreaks) ##You can provide your own intial values "par0" too
obj$par # the estimated parameters
obj$likelihood ##The likelihood at estimated parameters
obj$converge ##Whether the algorithm converged, "0" indicates convergence
obj$niter  ##Nubmer of iterations taken to converge

## Test for H0: g=0     
##Warning: if the algorithm does not converge, you need change values for par0 (initial values for alternative estimation) and par0_null (initial values for null estimation)
                 
obj <- gh_test(y,par0=NULL,par0_null=NULL,nbreaks=nbreaks,Test="g=0")
obj$STAT #the value of LRT statistics
obj$pvalue #The pvalue of the test


## Test for H0: h=0         
##Warning: if the algorithm does not converge, you need change values for par0 (initial values for alternative estimation) and par0_null (initial values for null estimation)
             
obj <- gh_test(y,par0=NULL,par0_null=NULL,nbreaks=nbreaks,Test="h=0")
obj$STAT #the value of LRT statistics
obj$pvalue #The pvalue of the test


## Test for H0: h=0, g=0    
                  
obj <- gh_test(y,par0=NULL,par0_null=NULL,nbreaks=nbreaks,Test="g=0,h=0")
obj$STAT #the value of LRT statistics
obj$pvalue #The pvalue of the test


######Regression model with g-h error y=x\beta+e###
X <- matrix(rnorm(n*5),n,5) ##Simulate design matrix
beta0 <- coef(lm(y~X-1))
res <- y-X%*%beta0
par0 <- Initials(res)
obj <- gh_likest_reg(y,X=X,beta0=beta0,par0=par0,nbreaks=nbreaks) ##You can provide your own intial values "par0" too