# HW 3 loaddat.R JLCY, 19 Dec 2013

############################################################ 

# clear variables load data set different variables

{

    number = Qnumber

    library(sm)

    source("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\skew.r")

    # This code is to model a lag-1 Modified Nonparametric K-nearest neighbor
    # model for a monthly streamflow time series.

    test = matrix(scan("C:\\Users\\jfklee90\\Desktop\\CVEN 6833\\HW 3\\LeesFerry-monflows-1906-2006.txt"), 
        ncol = 13, byrow = T)

    year = test[, 1]
    LFmonflow = test[, 2:13]
    LFmonflow = LFmonflow/10^6  # convert to MAF 

    # Obs May flow
    May = LFmonflow[, 5]

    # expand sequence with 100 values in between min & max
    xeval = seq(min(May) - 0.25 * sd(May), max(May) + 0.25 * sd(May), length = 100)

    nevals = length(xeval)

    nsim = 250  # number of simulations 

    nsim1 = nsim + 1

    # initialize the matrices that store the statistics the first row contains
    # the statistic of the historical data.

    armean = matrix(0, nsim, number)
    arstdev = matrix(0, nsim, number)
    arvar = matrix(0, nsim, number)
    arcor = matrix(0, nsim, number)
    arskw = matrix(0, nsim, number)
    armax = matrix(0, nsim, number)
    armin = matrix(0, nsim, number)

    if (number == 13) {

        armean = matrix(0, nsim1, number)
        arvar = matrix(0, nsim1, number)
        arsd = matrix(0, nsim1, number)
        arcor = matrix(0, nsim1, number)
        arskw = matrix(0, nsim1, number)
        armax = matrix(0, nsim1, number)
        armin = matrix(0, nsim1, number)
    }

    ## Points to evaluate the May pdf

    May = LFmonflow[, 5]

    # expand sequence with 100 values in between min & max
    xeval = seq(min(May) - 0.25 * sd(May), max(May) + 0.25 * sd(May), length = 100)

    nevals = length(xeval)

    nyrs = length(LFmonflow[, 1])  #number of years of data for each month
    nyrs1 = nyrs - 1

    ## Store the May simulations and array for the PDF
    maysim = 1:nyrs
    simpdf = matrix(0, nrow = nsim, ncol = nevals)

    Annual = apply(LFmonflow, 1, sum)
    annxeval = seq(min(Annual) - 0.25 * sd(Annual), max(Annual) + 0.25 * sd(Annual), 
        length = 100)

    annsim = 1:nyrs
    annsimpdf = matrix(0, nrow = nsim, ncol = nevals)

    sim_LS = 1:nsim  # longest surplus
    sim_MS = 1:nsim  # maximum surplus

    sim_LD = 1:nsim  # longest drought
    sim_MD = 1:nsim  # maximum drought

    # Scale the data
    fmean = apply(LFmonflow, 2, mean)  # monthly mean
    fsdev = apply(LFmonflow, 2, sd)  # monthly standard deviation

    LFmonflow.std = t((t(LFmonflow) - fmean)/fsdev)

}

## Error: 找不到物件 'Qnumber'