#Моделирование 

library("MASS")
library("matlab")
library("matrixcalc")
library("Matrix")



Fixed.Mod(2,2,2,100)
Random.Mod(2,2,2,100)
print(H.gen(3,3,3))
Fixed.Matrix.Mod(4,4,4,1000)
Fixed.Matrix.Diff.J.Mod(3,c(2,4,3),3,400)
Random.Matrix.Diff.J.Mod.Lambda(3,3,3,1000)

#Проверка фиксированного эффекта взаимодействия. Моделирование распределения pvalue для проверки гипотезы о равенстве 
#нулю эффекта взаимодействия
Fixed.Mod <- function(I,J,T,n){
sigma1 <- 3
sigma <- 2
SSE_a <- c()
SSE_c <- c()

mu <- 500
alpha <- c(1:I)
alpha[I] <- -sum(alpha[1:(I-1)])
beta <- c(1:T)
gamma <- matrix(c(1:(I*T)),ncol = T,nrow = I)
for(i in c(1:I)) gamma[I,i] <- -sum(gamma[c(1:(I-1)),i])
for(t in c(1:(T-1))) gamma[t,T] <- -sum(gamma[t,c(1:(T-1))])
gamma[I,T] <- -sum(gamma[I,-T])

for(p in c(1:n)){
  e <- array(rnorm(I*J*T,0,sigma),dim=c(I,J,T))
  e1 <- array(rnorm(I*J,0,sigma1),dim = c(I,J))
  x <- array(0,dim = c(I,J,T))
  for(i in c(1:I)){
    for(j in c(1:J)){
      for(t in c(1:T)){
        x[i,j,t] <- e[i,j,t] + e1[i,j] + alpha[i] + beta[t] + mu
      }
    }
  }
  
  SSE <- 0
  for(i in c(1:I)){
    for(j in c(1:J)){
      for(t in c(1:T)) SSE <- SSE + (x[i,j,t] - mean(x[i,j,]) - mean(x[i,,t]) + mean(x[i,,]))^2
    }
  }
  SSE_a <- c(SSE_a,SSE)
  
  SSС <- 0
  for(i in c(1:I)){
    for(j in c(1:J)){
      for(t in c(1:T)) SSС <- SSС + (mean(x[i,,t]) - mean(x[i,,]) - mean(x[,,t]) + mean(x[,,]))^2
    }
  }
  
  SSE_c <- c(SSE_c,SSС)  
}


F <- (SSE_c/((T-1)*(I-1)))/(SSE_a/((T-1)*(I*J-I)))

sample <- pf(F,(T-1)*(I-1),((T-1)*(I*J-I)))
q.e.chisq <- sample[order(sample)]

q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)

qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
abline(0, 1, col = 2)
}

Random.Mod <- function(I,J,T,n){
  sigma1 <- 3
  sigma <- 2
  sigma_a <- 5
  sigma_b <- 4
  sigma_c <- 2.5
  SSE_a <- 0
  SSE_c <- 0
  
  mu <- 500
  
  
  for(p in c(1:n)){
    alpha <- rnorm(I,0,sigma_a)
    beta <- rnorm(T,0,sigma_b)
    gamma <- matrix(rnorm(I*T,0,sigma_ab),nrow = I,ncol = T)
    e <- array(rnorm(I*J*T,0,sigma),dim=c(I,J,T))
    e1 <- array(rnorm(I*J,0,sigma1),dim = c(I,J))
    x <- array(0,dim = c(I,J,T))
    for(i in c(1:I)){
      for(j in c(1:J)){
        for(t in c(1:T)){
          x[i,j,t] <- e[i,j,t] + e1[i,j] + alpha[i] + beta[t] + mu #+ gamma[i,t]
        }
      }
    }
    
    SSE <- 0
    for(i in c(1:I)){
      for(j in c(1:J)){
        for(t in c(1:T)) SSE <- SSE + (x[i,j,t] - mean(x[i,j,]) - mean(x[i,,t]) + mean(x[i,,]))^2
      }
    }
    SSE_a <- c(SSE_a,SSE)
    
    SSС <- 0
    for(i in c(1:I)){
      for(j in c(1:J)){
        for(t in c(1:T)) SSС <- SSС + (mean(x[i,,t]) - mean(x[i,,]) - mean(x[,,t]) + mean(x[,,]))^2
      }
    }
    
    SSE_c <- c(SSE_c,SSС)  
  }
  
  SSE_a <- SSE_a[-1]
  SSE_c <- SSE_c[-1]
  
  F <- (SSE_c/((T-1)*(I-1)))/(SSE_a/((T-1)*(I*J-I)))
  
  sample <- pf(F,(T-1)*(I-1),((T-1)*(I*J-I)))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
  abline(0, 1, col = 2)
}

#Генерация матрицы плана
H.gen <- function(I,J,T){
  H <- matrix(0,nrow = I*J*T,ncol = I*(T-1))
  for(i in c(1:(I-1))){
    for(t in c(1:(T-1))){
      for(j in c(1:J)){
        H[(i-1)*J*T+j+(t-1)*J,1+(t-1)*I] <- 1
        H[(i-1)*J*T+j+(t-1)*J,2+(t-1)*I+(i-1)] <- 1
        H[(i-1)*J*T+j+(T-1)*J,1+(t-1)*I] <- -1
        H[(i-1)*J*T+j+(T-1)*J,2+(t-1)*I+(i-1)] <- -1
      }
    }
  }
  
  for(t in c(1:(T-1))){
    for(j in c(1:J)){
      H[(I-1)*J*T+(t-1)*J+j,(t-1)*I+1] <- 1
      for(p in (c(1:(I-1))))
        H[(I-1)*J*T+(t-1)*J+j,(t-1)*I+1+p] <- -1
    }
  }
  
  for(j in c(1:J)){
    for(t in c(1:(T-1))){
      H[I*J*T-J+j,I*(t-1)+1] <- -1
      for(i in c(1:(I-1))){
        H[I*J*T-J+j,I*(t-1)+1 + i] <- 1
      }
    }
  }

  H  
}  

#Проверка гипотезы \gamma = 0 для фиксированных эффектов в матричном виде
Fixed.Matrix.Mod <- function(I,J,T,n){
  sigma1 <- 3
  sigma <- 2
  sigma_a <- 5
  sigma_b <- 4
  sigma_ab <- 2.5
  SSE_a <- 0
  SSE_c <- 0
  H <- H.gen(I,J,T)
  
  mu <- 50
  
  Qe <- c()
  Qg <- c()
  
  alpha <- c(1:I)
  alpha[I] <- -sum(alpha[1:(I-1)])
  beta <- c(1:T)
  gamma <- matrix(c(1:(I*T)),ncol = T,nrow = I)
  for(i in c(1:I)) gamma[I,i] <- -sum(gamma[c(1:(I-1)),i])
  for(t in c(1:(T-1))) gamma[t,T] <- -sum(gamma[t,c(1:(T-1))])
  gamma[I,T] <- -sum(gamma[I,-T])
  
  
  
  for(p in c(1:n)){
    e <- array(rnorm(I*J*T,0,sigma),dim=c(I,J,T))
    e1 <- array(rnorm(I*J,0,sigma1),dim = c(I,J))
    x <- array(0,dim = c(I,J,T))
    for(i in c(1:I)){
      for(j in c(1:J)){
        for(t in c(1:T)){
          x[i,j,t] <- e[i,j,t] + e1[i,j] + alpha[i] + beta[t] + mu# + gamma[i,t]
        }
      }
    }
    
    Y <- c(1:(I*J*T))
    for(i in c(1:I)){
      for(t in c(1:T)){
        for(j in c(1:J)) {
          Y[j+J*(t-1)+T*J*(i-1)] <- x[i,j,t]-mean(x[i,j,])
        }
      }
    } 
    
    Qe <- c(Qe,t(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y)%*%(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y))
    
    Hb <- H[,c(0:(T-2))*I+1]
    Qg <- c(Qg,t(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y)%*%(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y))
  }
  
  F <- ((Qg-Qe)/((T-1)*(I-1)))/(Qe/(I*J*(T-1)-I*(T-1)))
  
  sample <- pf(F,(T-1)*(I-1),(I*J*(T-1)-I*(T-1)))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
  abline(0, 1, col = 2)
}

#Проверка гипотезы \sigma_\gamma = 0 для случайных эффектов в матричном виде
Random.Matrix.Mod <- function(I,J,T,n){
  sigma1 <- 3
  sigma <- 2
  sigma_a <- 5
  sigma_b <- 4
  sigma_ab <- 2.5
  SSE_a <- 0
  SSE_c <- 0
  H <- H.gen(I,J,T)
  
  mu <- 50
  
  Qe <- c()
  Qg <- c()
  
  for(p in c(1:n)){
    alpha <- rnorm(I,0,sigma_a)
    beta <- rnorm(T,0,sigma_b)
    gamma <- matrix(rnorm(I*T,0,sigma_ab),nrow = I,ncol = T)
    e <- array(rnorm(I*J*T,0,sigma),dim=c(I,J,T))
    e1 <- array(rnorm(I*J,0,sigma1),dim = c(I,J))
    x <- array(0,dim = c(I,J,T))
    for(i in c(1:I)){
      for(j in c(1:J)){
        for(t in c(1:T)){
          x[i,j,t] <- e[i,j,t] + e1[i,j] + alpha[i] + beta[t] + mu# + gamma[i,t]
        }
      }
    }
    
    Y <- c(1:(I*J*T))
    for(i in c(1:I)){
      for(t in c(1:T)){
        for(j in c(1:J)) {
          Y[j+J*(t-1)+T*J*(i-1)] <- x[i,j,t]-mean(x[i,j,])
        }
      }
    } 
    
    Qe <- c(Qe,t(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y)%*%(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y))
    
    Hb <- H[,c(0:(T-2))*I+1]
    Qg <- c(Qg,t(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y)%*%(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y))
  }
  
  F <- ((Qg-Qe)/((T-1)*(I-1)))/(Qe/(I*J*(T-1)-I*(T-1)))
  
  sample <- pf(F,(T-1)*(I-1),(I*J*(T-1)-I*(T-1)))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
  abline(0, 1, col = 2)
}

#Проверка гипотезы \sigma_\gamma = 0 для случайных эффектов в матричном виде для разных J
Fixed.Matrix.Diff.J.Mod <- function(I,J,T,n){
  
  Hj.gen <- function(H,I,J,T,Jm){
    answ <- c()
    for(i in c(1:I))
      for(t in c(1:T)){
        for(j in c(1:J)){
          curr <- j+J*(t-1)+J*T*(i-1)
          if(mod(curr-1,J)<=Jm[i]-1) answ <- c(answ,curr)
        }
      }
    answ
  }
  
  Jl <- length(J)
  sigma1 <- 3
  sigma <- 2
  sigma_a <- 5
  sigma_b <- 4
  sigma_ab <- 2.5
  SSE_a <- 0
  SSE_c <- 0
  J_max <- max(J)
  H <- H.gen(3,J_max,3)
  J_count <- c(2,4,3)
  J <- Hj.gen(H,3,4,3,J_count)
  J_n <- length(J)
  H <- H[J,]
  mu <- 500
  Qe <- c()
  Qg <- c()
  
  alpha <- c(1:I)
  alpha[I] <- -sum(alpha[1:(I-1)])
  beta <- c(1:T)
  gamma <- matrix(c(1:(I*T)),ncol = T,nrow = I)
  for(i in c(1:I)) gamma[I,i] <- -sum(gamma[c(1:(I-1)),i])
  for(t in c(1:(T-1))) gamma[t,T] <- -sum(gamma[t,c(1:(T-1))])
  gamma[I,T] <- -sum(gamma[I,-T])
  
  for(p in c(1:n)){
    e <- array(rnorm(I*J_max*T,0,sigma),dim=c(I,J_max,T))
    e1 <- array(rnorm(I*J_max,0,sigma1),dim = c(I,J_max))
    x <- array(0,dim = c(I,J_max,T))
    for(i in c(1:I)){
      for(j in c(1:J_count[i])){
        for(t in c(1:T)){
          x[i,j,t] <- e[i,j,t] + e1[i,j] + alpha[i] + beta[t] + mu# + gamma[i,t]
        }
      }
    }
    
    Y <- c(1:(I*J_max*T))
    for(i in c(1:I)){
      for(t in c(1:T)){
        for(j in c(1:J_max)) {
          Y[j+J_max*(t-1)+T*J_max*(i-1)] <- x[i,j,t]-mean(x[i,j,])
        }
      }
    } 
    
    Y <- Y[J]
    
    Qe <- c(Qe,t(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y)%*%(Y - H%*%ginv(t(H)%*%H)%*%t(H)%*%Y))
    
    Hb <- H[,c(0:(T-2))*I+1]
    Qg <- c(Qg,t(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y)%*%(Y - Hb%*%ginv(t(Hb)%*%Hb)%*%t(Hb)%*%Y))
  }
  
  F <- ((Qg-Qe)/((T-1)*(I-1)))/(Qe/(J_n-sum(J_count)-I*(T-1)))
  
  sample <- pf(F,(T-1)*(I-1),(J_n-sum(J_count)-I*(T-1)))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
  abline(0, 1, col = 2)
}

#Проверка гипотезы \sigma_\gamma = 0 для случайных эффектов в матричном виде с недиагональной матрицей ковариаций ошибок
Random.Matrix.Diff.J.Mod.Lambda <- function(I,J,T,n){
  sigma1 <- 3
  sigma <- 1
  sigma_a <- 5
  sigma_b <- 4
  sigma_ab <- 2.5
  SSE_a <- 0
  SSE_c <- 0
  H <- H.gen(I,J,T)
  
  mu <- 1000
  
  Qe <- c()
  Qg <- c()
  
  for(p in c(1:n)){
    N <- I*J*T
    Lambda <- matrix(runif(N*N,1,4),nrow=N,ncol=N)
    for(i in c(1:N))
      for(j in c(i:N)){
        Lambda[j,i] <- Lambda[i,j]  
      }
    alpha <- rnorm(I,0,sigma_a)
    beta <- rnorm(T,0,sigma_b)
    gamma <- matrix(rnorm(I*T,0,sigma_ab),nrow = I,ncol = T)
    e <- Lambda%*%array(rnorm(I*J*T,0,sigma),dim=c(I,J,T))
    e1 <- array(rnorm(I*J,0,sigma1),dim = c(I,J))
    
    
    Y <- c(1:(I*J*T))
    for(i in c(1:I)){
      for(t in c(1:T)){
        for(j in c(1:J)) {
          Y[j+J*(t-1)+T*J*(i-1)] <- e[j+J*(t-1)+T*J*(i-1)]+beta[t]#+gamma[i,t]
        }
      }
    } 
    
    
    Lambda <- Lambda%*%t(Lambda)
    Qe <- c(Qe,t(Y - H%*%ginv(t(H)%*%ginv(Lambda)%*%H)%*%t(H)%*%ginv(Lambda)%*%Y)%*%ginv(Lambda)%*%(Y - H%*%ginv(t(H)%*%ginv(Lambda)%*%H)%*%t(H)%*%ginv(Lambda)%*%Y))
    #Qe <- c(Qe,t(e)%*%ginv(Lambda)%*%e)
    
    Hb <- H[,c(0:(T-2))*I+1]
    Qg <- c(Qg,t(Y - Hb%*%ginv(t(Hb)%*%ginv(Lambda)%*%Hb)%*%t(Hb)%*%ginv(Lambda)%*%Y)%*%ginv(Lambda)%*%(Y - Hb%*%ginv(t(Hb)%*%ginv(Lambda)%*%Hb)%*%t(Hb)%*%ginv(Lambda)%*%Y))
  }
  mean(Qg-Qe)
  
  sample <- pchisq(Qe,I*J*T-I*(T-1))
  sample <- pchisq(Qg,I*J*T - T + 1)
  #sample <- pchisq(Qg-Qe,(I-1)*(T-1))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Empirical quantiles", ylab = "Theoretical quantiles")
  abline(0, 1, col = 2)
  
  
  F <- ((Qg-Qe)/((T-1)*(I-1)))/(Qe/(I*J*T-I*(T-1)))
  
  sample <- pf(F,(T-1)*(I-1),(I*J*T-I*(T-1)))
  q.e.chisq <- sample[order(sample)]
  
  q.unif <- sapply(0:(n-1)/n,FUN = qunif,min=0,max=1)
  
  qqplot(q.e.chisq,q.unif,xlab = "Эмпирические квантили", ylab = "Теоретические квантили")
  abline(0, 1, col = 2)
}