#NACITANIE DAT
Yt <- c(6.1, 7.3, 9.6, 10.2, 10.1, 11.3, 12.2, 12.5, 13.2)
Pt <- c(103, 102, 100, 94, 98, 97, 98, 97, 96)
It <- c(110, 114, 130, 135, 141, 152, 160, 165, 170)

#príprava matíc
mat=cbind(Yt,Pt,It)
model<-(Yt ~ Pt + It)
regresia<-lm(model)

#Goldfield-Quandt Test
gqtest(regresia, order.by = Pt, fraction=1, data = mat)

#Breusch-Pagan Test
bptest(regresia)

#Bartlett test
data<-list(Yt,Pt,It)
bartlett.test(data)
Yt1<-1

for (i in 1:9)
{
  Yt1[i]<-(Yt[i]/It[i])
}

Xt1<-1

for (i in 1:9)
{
  Xt1[i]<-(Pt[i]/It[i])
}
mat1 = cbind(Yt1, Xt1, It)
model1<-(Yt1 ~ Xt1 + It)
regresia1<-lm(model1)
gqtest(regresia1, order.by = It, fraction=2, data = mat1)

#VLASTNA FUNKCIA
Het = function (mat, x, m,rem)
{
  nr = nrow(mat)
  nc = ncol(mat)
  mat1 = mat
  matt=mat
  mat1 = mat1[order(mat1[,x+1]),]  
  if (nr%%2==1){
    if (m%%2==1){
      k=m
    }  
    else{
      k=m+1
    }
  } 
  else {
    if (m%%2==1){
      k=m+1
    }  
    else{
      k=m
    }
  }
  l=0
  l=(nr-k)/2
  mt1 = mat1[1:l,1:nc] 
  cat(mt1,"\n")
  mt2 = mat1[(l+k+1):nr,1:nc]
  cat(mt2,"\n")
  y1=mt1[1:l,1]
  cat(y1,"\n")
  y2=mt2[1:l,1]
  cat(y2,"\n")
  reg = function (Yt, mm)
  {
    regr = lm (Yt~mm)
    nr = nrow(mm)
    nc = ncol(mm)
    y = c()
    b = c()
    for (i in 1:(nc+1))
    {
      b[i]=regr$coef[i]
    }    
    for (i in 1:nr)
    {
      ypom = 0;
      for (j in 2:(nc+1))
      {
        ypom = ypom + b[j]*mm[i,(j-1)]
      }
      y[i]=b[1]+ypom
    }
    y   
  }
  yt1 = reg(y1,mt1[1:l,2:nc])
  yt2 = reg(y2,mt2[1:l,2:nc])
  sum1 = 0
  sum2 = 0
  et1=c()
  et2=c()
  for (i in 1:l){
    et1[i]=y1[i]-yt1[i]
    et1[i]=et1[i]*et1[i]
    sum1 =sum1 + et1[i]
  }
  for (i in 1:l){
    et2[i]=y2[i]-yt2[i]
    et2[i]=et2[i]*et2[i]
    sum2 = sum2 + et2[i]
  }
  if (sum1>sum2){
    Fv=sum1/sum2
  }
  else {
    Fv=sum2/sum1
  }
  sv = l-(nc-1)-1
  Ftab = qf(0.95,sv,sv)
  if (Fv>Ftab){
    H1='true'
    cat(
      " =======================================","\n",
      " Results:","\n\n",
      "---------------------------------------","\n",
      " F stats. calculated: ",Fv,"\n",
      " F stats. tabulated: ",Ftab,"\n",
      " Presence of heteroscedasticity","\n",
      "=======================================","\n\n"
    )
  }
  else {
    H1='false'
    cat(
      " =======================================","\n",
      " Results:","\n\n",
      "---------------------------------------","\n",
      " F stats. calculated: ",Fv,"\n",
      " F stats. tabulated: ",Ftab,"\n",
      " The absence of heteroskedasticity","\n",
      "=======================================","\n\n"
    )
  }
  
  if (H1=='true' & rem=='true'){
    cat(
      " +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+","\n",
      " Removing heteroscedasticity","\n",
      "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+","\n\n"
    )
    xx=matt[1:nr,x+1]
    for (i in 1:nc){
      if (i==x+1){
        for (j in 1:nr){
          matt[j,i]=1/xx[j]
        }
      }
      else{
        for (j in 1:nr){
          matt[j,i]=matt[j,i]/xx[j]
        }
      }
    }   
    Het(matt,1,1,'false')  
  }
  if (rem=='true'){
    cat(
      " =======================================","\n",
      " New values:","\n\n",
      "---------------------------------------","\n",
      "yt=c("
    )
    for (j in 1:nr){
      if(j<nr){
        cat(
          matt[j,1],","
        )
      }
      else{
        cat(
          matt[j,1]
        )
      }
    }
    cat(
      ")\n"
    )
    for (i in 2:nc){
      cat(
        " x",i-1,"=",matt[1:nr,i],"\n"
      )
    }     
    cat(
      " =======================================","\n\n"
    )
  }
}
Het(cbind(Yt,Pt,It),1,1,'true')